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A.  SCIENTIFIC  OBJECTIVES: 

Conventional  radar  signal  processing  is  based  on  two  simplified  assumptions  about  target 
scattering:  (i)  that  the  target  is  a  rigid  body;  (ii)  that  the  target  can  be  modeled  as  a 
collection  of  independent  point  scatterers  without  any  multiple  scattermg  effects. 
However,  real  radar  data  can  deviate  significantly  from  these  two  simplified  assumptions. 
First,  real-world  targets  are  often  observed  by  radar  sensors  under  dynamic  conditions 
where  non-rigid  body  motions  can  exist.  These  non-rigid  body  motions  give  rise  to 
“microDoppler”  phenomena,  which  have  been  observed  in  a  number  of  SAR  and  ISAR 
sensors.  Examples  of  microDoppler  phenomena  include  returns  from  moving 
components  on  the  target  such  as  scanning  antennas  or  rotating  wheels,  as  well  as  those 
from  flexing  and  vibration  of  the  target  frame.  Second,  strong  multiple  scattering  physics 
are  often  encountered  in  inlets  and  cavity  structures  on  the  target.  For  instance,  the  most 
prominent  feature  on  an  air  target  is  often  the  range-delayed  return  from  the  jet  inlet  duct. 
Significant  modeling  work  has  been  carried  out  by  the  computational  electromagnetics 
community  to  characterize  the  complex  scattering  from  inlet  structures,  yet  little  effort 
has  been  placed  on  utilizing  the  results  to  develop  better  imaging  algorithms  to  map  the 
inlet  interior.  The  objectives  of  this  research  program  are:  (i)  to  gam  in-depth 
understanding  of  these  higher-order  phenomena  through  simulation  and  measurement,  (ii) 
to  develop  physics-based  models  and  the  associated  signal  processing  strategies  to  extract 
the  resulting  radar  features,  and  (iii)  to  exploit  and  utilize  these  additional  features  to 
enhance  the  performance  of  automatic  target  recognition  (ATR)  algorithms. 


B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS: 

During  the  first  year  of  this  research  program,  three  areas  of  investigation  have  been 
initiated.  First,  we  have  designed,  built  and  tested  a  low-cost  radar  sensor  for 
microDoppler  data  collection.  The  system  costs  less  than  $3,000  and  is  tunable  between 
4  to  10  GHz.  Preliminary  data  from  a  moving  human  have  been  collected  and  processed. 
Second,  we  have  begun  to  develop  simulation  tools  to  aid  in  the  interpretation  and 
understanding  of  microDoppler  and  multiple  scattering  phenomena  in  complex  targets. 
In  particular,  a  multi-platen  z-buffer  algorithm  was  investigated  for  fast  ray  tracing  and 
the  results  show  promise  as  a  real-time  simulation  tool.  Finally,  we  have  investigated  the 
application  of  genetic  algorithms  and  time-reversal  imaging  techniques  to  improve  target- 
to-clutter  performance  and  to  achieve  better  feature  extraction.  Our  detailed  progress 
along  these  three  lines  is  described  below. 

MicroDoppler  Radar  and  MicroDoppler  Data  Collection.  We  have  developed  a 
wideband  radar  testbed  that  is  fully  tunable  between  4  to  10  GHz  (see  Fig.  1).  The 
purpose  of  this  hardware  effort  is  to  provide  data  collection  capability  that  will  allow  us 
to  study  the  microDoppler  phenomenology  in  great  detail.  The  radar  is  a  homodyne 
system  with  full  I/Q  detection.  Data  acquisition  is  accomplished  using  a  National 
Instruments  DAQ  card  and  the  results  can  be  displayed  on  a  laptop  computer.  Careful 
calibration  of  the  radar  was  carried  out  using  tuning  forks  and  an  audio  speaker. 
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Preliminary  data  collection  of  walking  humans  was  carried  out  in  an  indoor,  high- 
clutter  environment  (see  Fig  2).  The  microDoppler  phenomenon  from  human  gait  was 
previously  investigated  by  Geisheimer  et  al.  [1,  2].  It  was  shown  that  a  number  of 
interesting  microDoppler  features  from  a  walking  person  can  be  observed  using  the 
spectrogram.  A  similar  analysis  was  also  carried  out  by  processing  the  data  collected 
from  Navy’s  APY-6  radar  of  a  person  walking  [3].  However,  much  work  is  needed  to 
gain  a  more  complete  phenomenological  understanding  of  the  detailed  features  associated 
with  human  movements  in  a  wide  variety  of  scenarios. 
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Fig.  2.  MicroDoppler  features  observed  in  the  measured  data  of  a  walking  person 
collected  using  the  radar  shown  in  Fig.  1. 

We  have  carried  out  some  preliminary  feature  extraction  effort  using  the  joint 
time-frequency  (JTF)  processing  tools  that  we  have  developed  previously  under  ONR 
sponsorship.  Fig.  3  shows  the  use  of  chirplet  bases  to  achieve  signal  separation  of  the 
human  gait  data.  The  different  (dwell  time)-(Doppler  frequency)  behaviors  from  various 
body  movements  were  exploited  to  separate  out  the  individual  microDoppler  features. 

In  the  coming  year,  we  plan  to  carry  out  extensive  data  collection  over  a  wide 
variety  of  parameter  spaces,  include  various  types  of  movements  (walking,  running, 
crawling,  and  involuntary  motions  such  as  respiration),  different  environments  (indoor, 


Fig.  3.  Feature  extraction  of  the  measured  walking  data  using  joint  time-frequency  signal 
processing.  The  raw  data  are  parameterized  using  the  adaptive  chirplet 
decomposition,  and  the  different  features  are  extracted  based  on  the  chirplet 
parameters. 

outdoor  urban,  through-wall,  occluded  by  foliage,  underground  facilities),  over  a  wide 
range  of  frequencies,  and  for  different  human  subjects.  We  will  investigate  the  use  of 
improved  basis  sets  derived  from  the  actual  physics  of  human  movements  for  feature 
extraction.  As  this  research  matures,  we  expect  to  build  up  a  comprehensive  dictionary 
of  microDoppler  features,  from  which  feature  extraction  and  automatic  classification 
algorithms  can  be  successfully  developed. 

Ray  Tracing  Simulation.  We  have  initiated  research  to  develop  efficient  simulation 
tools  to  aid  in  the  interpretation  and  understanding  of  microDoppler  and  multiple 
scattering  phenomena  in  complex  targets.  Our  simulation  methodology  is  based  on  the 
shooting  and  bouncing  ray  (SBR)  technique.  The  simulation  of  radar  returns  from  non- 
rigid  targets  is  extremely  computationally  intensive  based  on  SBR.  Therefore,  we  have 
devoted  our  efforts  to  explore  ways  to  speed  up  the  ray  tracing  time. 

The  Binary  Space  Partition  (BSP)  tree  algorithm  [4]  is  the  most  standard  ray 
tracer  in  use.  In  the  algorithm,  a  BSP  tree  is  first  built  based  on  the  facet  model  of  the 


target  by  recursively  cutting  the  bounding  box  of  the  object  along  a  spatial  plane.  Ray 
tracing  is  then  performed  by  traversing  the  BSP  tree.  The  BSP-tree  based  ray  tracer  is 
considered  the  fastest  among  all  of  the  spatial  subdivision  approaches.  Recently,  the 
multiplaten  z-buffer  ray-tracing  algorithm  was  proposed  by  Hu  et  al  [5,  6]  as  an 
alternative  to  the  traditional  BSP  tree  algorithm.  In  the  multiplaten  z-buffer  (MPZ) 
approach,  a  multi-layered  z-buffer  is  first  generated  from  the  scan  conversion  process 
(Fig.  4(a)).  Instead  of  just  storing  the  z-coordinates  of  the  visible  pixels  as  in  the 
traditional  z-buffer  process,  multiple  z-buffers  are  created  to  store  the  z-coordinates  of  all 
of  the  facets  within  each  pixel  during  the  scan  conversion.  During  the  ray  trace,  a  ray  is 
tracked  by  moving  along  the  ray  direction  pixel-by-pixel.  Within  every  pixel,  the  z-depth 
of  the  ray  is  compared  to  all  of  the  z-buffer  values  for  that  pixel  to  check  for  possible 
intersections  (Fig.  4(b)).  Once  an  intersection  is  found,  the  hit  point  and  the  reflection 
direction  can  be  calculated,  and  the  tracing  process  is  then  iterated  until  the  ray  departs 
from  the  bounding  box.  We  have  evaluated  the  computation  time  performance  of  the 
MPZ  ray  tracer  against  that  of  the  BSP  tree-based  algorithm.  Results  for  a  wide  range  of 
targets  have  been  tested  to  determine  the  computational  as  well  as  memory  complexity  as 
functions  of  the  number  of  facets  and  the  complexity  of  the  target.  We  found  that, 
contrary  to  the  BSP  algorithm,  the  complexity  of  the  MPZ  is  independent  of  the  number 
of  facets  comprising  the  target  (see  Fig.  5).  However,  the  computation  time  is  dependent 
on  the  number  of  pixels  traversed  by  a  ray. 


Multi-layered  Z-buffer 


Fig.  4.  (a)  Multi-layered  z-buffer.  (b)  Ray  tracing  using  the  multilayered  z-buffer. 
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Fig.  5.  Computational  time  of  the  various  ray-tracing  algorithms  versus  the  number  of 
facets  comprising  a  test  cavity,  (a)  Binary  Space  Partition  (BSP)  tree,  (b)  Multi¬ 
platen  z-buffer  (MPZ).  (c)  Multi-aspect  MPZ  (MMPZ)  with  a  z-buffer  every  20°. 
(d)  Multi-aspect  MPZ  with  a  z-buffer  every  1°. 


To  further  speed  up  the  performance  of  MPZ  algorithm,  we  are  investigating  a 
multi-aspect  MPZ  approach.  The  approach  is  motivated  by  the  fact  that  the  number  of 
pixels  a  ray  traverses  between  bounces  can  be  reduced  dramatically  by  decreasing  the 
angle  between  the  ray  direction  and  the  z-buffer  direction.  In  the  algorithm,  multiple 
multi-layered  z-buffers  are  first  generated  from  the  scan  conversion  process  along  many 
aspect  angles.  The  maximum  number  of  multi-layered  z-buffers  is  limited  only  by  the 
available  memory  resource.  The  more  aspect  angles  that  be  stored,  the  less  pixels  a  ray 
traverses  in  one  bounce,  and  the  better  the  time  performance.  During  the  ray  trace,  the 
multi-layered  z-buffer  structure  that  has  the  closest  aspect  to  the  ray  direction  is  selected 
to  carry  out  the  ray  tracing.  A  ray  is  then  tracked  by  moving  along  the  ray  direction  inside 
this  MPZ  structure  pixel-by-pixel  to  check  for  possible  intersections.  Once  an  intersection 
is  found,  the  hit  point  and  the  reflection  direction  are  calculated.  Based  on  the  new  ray 
direction,  a  new  MPZ  is  chosen,  and  the  tracing  process  is  iterated  until  the  ray  departs 


from  the  bounding  box.  As  can  be  seen  in  Fig.  5,  a  dramatic  improvement  of  performance 
against  that  of  the  single-aspect  MPZ  can  be  achieved.  We  are  currently  investigating  the 
possibility  of  implementing  this  algorithm  as  a  real-time  tool  to  simulate  scattering  from 

complex  targets. 

Exploratory  Research  on  Inverse  Scattering  and  Time-Reversal  Imaging.  We  have 
carried  out  some  exploratory  research  into  the  topics  of  inverse  scattering  and  time- 
reversal  imaging.  It  is  well  recognized  that  traditional  imaging  algorithms  suffer  from 
resolution  limitation  and  image  artifacts  due  to  multiple  scattering  phenomena. 
Rigorously  solving  the  electromagnetic  inverse  scattering  problem,  on  the  other  hand,  is 
much  more  challenging.  In  our  previous  research  on  inverse  scattering,  we  applied  a 
genetic  algorithm  (GA)  together  with  a  computational  electromagnetics  solver  to  attack 
the  two-dimensional  inverse  scattering  problem  [7].  The  inversion  was  cast  into  an 
optimization  problem  whereby  the  difference  between  the  measured  fields  and  the 
computed  fields  from  a  forward  electromagnetic  solver  was  minimized.  We  found  that 
while  GA  was  well  suited  in  searching  for  the  global  optimum,  it  suffered  from  slow 
convergence.  Since  the  evolutionary  process  for  the  standard  GA  to  reach  a  cost 
minimum  is  in  general  very  slow  in  comparison  to  a  local  search  algorithm,  a  natural 
improvement  to  speed  up  the  simple  GA  is  to  hybridize  the  simple  GA  with  a  local 
search.  While  this  hybrid  GA  (HGA)  shows  improvements  in  performance,  it  also  leads 
to  some  inefficiency.  As  the  parent  selection  scheme  of  GA  gives  priority  to  the  best 
members,  it  usually  leads  to  a  population  that  is  highly  clustered  around  the  local  minima. 
This  clustering  is  necessary  for  the  simple  GA  to  evolve  closer  to  the  exact  minimum.  For 
HGA,  however,  since  the  local  minima  have  been  completely  explored  by  the  local 
search,  such  clustering  will  lead  to  the  re-exploration  of  these  regions,  which  is  quite 
wasteful.  In  our  research,  we  have  proposed  a  technique  combining  HGA  with  the  tabu 
list  concept  to  increase  the  efficiency  of  the  HGA.  The  idea  was  motivated  by  the  Tabu 
Search  (TS)  algorithm  used  in  combinatorial  problems  [8,  9].  The  tabu  list  is  adopted  to 
exclude  those  regions  in  the  parameter  space  that  have  already  been  explored  by  the  local 
search.  In  this  manner,  there  will  be  no  revisiting  of  the  explored  regions  and  the  GA 
population  can  be  spread  out  to  explore  new  regions,  thus  improving  the  search 


efficiency.  We  have  applied  this  algorithm  to  the  electromagnetic  inverse  problem  of 
shape  reconstruction  of  metallic  cavity  structures  containing  strong  multiple  scattering 
effects.  The  algorithm  has  been  applied  to  reconstruct  the  shape  of  a  metallic  cavity 
based  on  the  measured  Ipswich  data  [10].  Fig.  6  shows  that  the  inversion  results  from  the 
HGA-Tabu  converge  faster  and  at  higher  success  rate  than  those  of  the  simple  GA  and 
hybrid  GA.  The  algorithm  could  potentially  be  useful  in  other  optimization/inverse 

problems. 
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Fie  6.  (a)  Convergence  comparison  for  the  shape  inversion  of  a  metallic  cavity  for 
b  random  search,  simple  GA,  hybrid  GA  (HGA)  and  HGA-Tabu.  The  HGA-Tabu 
shows  the  best  convergence,  (b)  Typical  inverted  shape  based  on  measured  data 
by  simple  GA.  (c)  Typical  inverted  shape  by  HGA-Tabu. 


Finally,  we  have  initiated  some  exploratory  research  into  the  concept  of  time 
reversal  imaging.  In  the  past  decade,  time  reversal  methods  have  been  extensively 
investigated  in  the  acoustics  and  mathematics  communities  [11-14].  The  original 
motivation  of  time  reversal  was  to  achieve  lithotripsy  by  focusing  the  acoustic  field  from 
a  near-field,  physical  array  on  a  prescribed  focal  point  through  a  complex  propagation 
medium.  This  was  accomplished  by  applying  the  time  reversal  operation  (or  phase 
conjugation  in  frequency)  to  a  pilot  signal  sent  out  by  the  array.  It  was  demonstrated 
theoretically  and  experimentally  that  the  fields  from  the  array  could  be  focused  at  a 
prescribed  point  through  a  complex  propagation  medium.  Thus,  it  is  potentially  an 
attractive  means  of  overcoming  medium  distortion  and  clutter  in  radar  imaging  of 
occluded  targets.  However,  there  are  major  gaps  between  what  has  been  achieved  in 
acoustics  and  what  needs  to  be  developed  for  radar  sensors  to  achieve  high-resolution 
radar  imaging.  First,  it  has  not  been  demonstrated  to  date  that  being  able  to  focus  on  a 
point  scatterer  can  translate  directly  into  good  imaging  performance  without  additional 
information  about  the  medium.  Nor  is  it  clear  that  any  proposed  algorithms  can  perform 
effectively  in  realistic  radar  clutter  such  as  foliage  or  walls.  Second,  existing  time 
reversal  methods  require  the  collection  of  array  calibration  data  in  the  bistatic  scenario 
involving  a  real  physical  array.  Realistic  radar  imaging  sensors,  on  the  other  hand, 
operate  in  the  monostatic,  synthetic  aperture  mode.  Finally,  contrary  to  the  scalar  wave 
case  considered  in  acoustics,  electromagnetic  polarization  is  a  new  dimension  that  should 
be  exploited  in  time  reversal  imaging  for  clutter  rejection. 

In  our  preliminary  research,  we  address  the  case  when  the  medium  is  unknown, 
yet  a  set  of  “calibration  targets”  with  known  positions  is  located  in  the  vicinity  of  the 
target  to  be  imaged.  This  means  that  once  we  can  focus  the  wave  on  the  calibration  target 
using  time  reversal,  it  is  possible  to  extract  useful  information  about  the  propagation 
distortion  caused  by  the  clutter.  Fig.  7(a)  shows  the  scenario  considered.  We  use  point 
scatterer  testing  and  assume  that  the  medium  distortion  is  described  by  a  random  phase 
function  <2>that  is  uniformly  distributed  between  [0,2 n]  and  uncorrelated  between  paths. 
Fig.  7(b)  shows  the  results  of  applying  conventional  free-space  imaging  to  the  data.  The 
smearing  in  the  cross  range  is  severe.  Next,  we  apply  a  subspace-based  time  reversal 
algorithm  [12]  to  the  data  and  focus  the  wave  onto  a  calibration  point  (x0,y0)  •  After  the 
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(c)  Applying  time  reversal  imaging  with  a  (d)  Reference  clutter-free  image, 

known  calibration  point. 


Fig.  7.  Results  of  time  reversal  imaging  by  focusing  on  a  calibration  point. 


time  reversal  process,  the  synthesized  receive  fields  on  the  array  elements  are: 


E  ( kx > ^ 


-j(kxx0+kyyo)e-j<P 


(1) 


where  kx  and  ky  are  related  to  the  array  element  position  and  radar  frequency.  Since 
(x0,y0)  are  known,  we  can  now  remove  the  propagation  effect  and  construct  the  image. 

The  result  is  shown  in  Fig.  7(c).  Fig.  7(d)  shows  the  reference  image  in  the  absence  of 
any  clutter.  As  we  can  see,  the  time  reversal  processing  with  calibration  nearly  recovers 
the  clutter-free  image.  The  minor  “ghosts”  are  believed  to  be  due  to  the  assumption  that 
the  eigenvector  chosen  corresponds  only  to  the  calibration  point.  This  deficiency  can  be 
circumvented  by  carrying  out  a  MUSIC-like  procedure  using  null-space  processing  [15]. 
There  are  many  open  questions  still  yet  to  be  addressed  theoretically,  and  much  more 
research  is  needed  to  make  this  concept  a  practical  radar  algorithm. 


C.  FOLLOW-UP  STATEMENT: 

During  the  coming  year,  our  research  will  be  devoted  to  three  fronts.  First,  we 
will  utilize  our  microDoppler  radar  testbed  to  carry  out  extensive  data  collection. 
Various  scenarios  will  be  considered,  including  types  of  movements,  different 
environments,  over  a  wide  range  of  frequencies,  and  for  different  targets.  Second,  we 
will  apply  simulation  and  signal  processing  tools  for  the  extraction  and  interpretation  of 
scattering  phenomenology  in  the  measured  data.  Third,  we  will  leverage  upon  the 
measured  database  to  develop  classification  algorithms  to  distinguish  different  classes  of 
targets  to  achieve  automatic  recognition.  The  potential  impacts  of  this  basic  research 
program  on  naval  radar  surveillance  capabilities  are  twofold.  First,  through  simulation, 
measurement  and  signal  processing,  this  research  will  result  in  an  in-depth  understanding 
of  the  radar  features  due  to  target  micro-dynamics  and  multiple  scattering  physics. 
Second,  our  research  will  lead  to  novel  exploitation  of  target  features  that  can 
significantly  improve  the  non-cooperative  target  recognition  (NCTR)  capabilities  for 
existing  and  future  naval  radar  sensors. 
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Abstract  We  present  an  algorithm  to  detect  the  presence  of  3D  target  motion  from  ISAR  data.  Based  on  the  3D 
point  scatterer  model,  we  first  examine  the  effect  of  3D  motion  on  ISAR  imaging.  It  is  shown  that  existing  motion 
compensation  algorithms  cannot  properly  focus  targets  exhibiting  3D  motion  during  the  imaging 
algorithm  is  then  derived  to  blindly  detect  the  degree  of  3D  target  mot.on  from  raw  radar  data.  It  M  on 
measuring  the  linearity  of  phases  between  two  or  more  point  scattered  on  the  target.  The  phase  esttmat  on  is 
implemented  using  the  adaptive  joint  time-frequency  technique.  Examples  are  provided  to  demonstra  e 
effectiveness  of  the  3D  motion  detection  algorithm  with  both  simulation  and  real  ISAR  data.  The  detection  results 
are  corroborated  with  the  truth  motion  data  from  on-board  motion  sensors  and  correlated  with  the  resulting 

ISAR  images. 
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1.  Introduction 

High-resolution  inverse  synthetic  aperture  radar  (ISAR)  imaging  is  regarded  as  an 
effective  tool  in  automatic  target  recognition  [1],  [2].  Ideally,  the  desired  target  motion 
is  uniform  rotation  without  translational  motion,  under  which  a  two-dimensional  (2D) 
Fourier  transform  brings  the  radar  data  in  the  (frequencyMdwell  time)  domain  into  the 
(range)-(Doppler  frequency)  domain.  Otherwise,  motion  compensation  is  needed  as  an 

intermediate  step  to  form  a  focused  ISAR  image.  .  t 

Since  target  motion  can  always  be  decomposed  into  translational  motion  and  rotational 
motion,  a  typical  motion  compensation  algorithm  consists  of  two  steps.  First,  a  point  on 
the  target  is  focused  through  translational  motion  compensation.  When  there  is  non- 
uniform  rotational  motion,  other  points  on  the  target  are  not  necessarily  focused. 
Rotational  motion  compensation  is  then  applied  to  focus  these  other  points.  Existing 
motion  compensation  algorithms  usually  assume  that  the  rotational  motion  of  a  target  is 
confined  to  a  2D  plane  during  the  dwell  duration  [1],  [2],  [3],  [4],  [5],  [6],  [7].  We  shall 
use  the  term  2D  motion  to  refer  to  target  rotational  motion  of  this  type.  Under  the  2D 
motion  assumption,  rotational  compensation  of  a  second  point  on  the  target  will  focus 
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the  whole  target.  When  there  is  3D  motion ,  i.e.,  when  the  rotational  motion  is  not 
confined  to  a  2D  plane,  rotational  compensation  of  a  second  point  cannot  focus  the 
whole  target. 

Recently,  several  independent  research  groups  have  reported  that,  for  aircraft  undergoing 
fast  maneuvers  or  ships  on  rough  seas,  the  motion  of  a  target  may  be  highly  chaotic  and 
does  not  always  obey  the  2D  motion  model  [8],  [9],  [10].  As  a  result,  the  image  formed 
using  the  standard  motion  compensation  algorithms  is  blurred.  In  [8]  and  [9],  the  effect  of 
3D  motion  on  ISAR  imaging  is  discussed.  However,  target  motions  are  assumed  to  be 
known  from  other  auxiliary  sensor  data  that  are  usually  not  accessible  in  real  operational 
environment  In  [10],  the  imaging  interval  is  adaptively  chosen  based  on  the  resolved 
target  feature  in  the  radar  image  to  overcome  the  3D  motion  issue.  It  requires  sound 
knowledge  of  the  target  under  consideration,  which  is  often  not  known  to  the  end  users  of 
ISAR  data. 

The  objective  of  this  paper  is  to  develop  an  algorithm  to  detect  the  presence  of  3D 
motions  during  the  imaging  interval  from  ISAR  data.  Based  on  the  3D  point  scatterer 
model,  we  first  examine  the  effect  of  3D  motion  on  existing  imaging  algorithms.  We 
then  develop  an  algorithm  to  blindly  detect  the  existence  of  3D  motion.  For  this 
purpose,  only  the  estimation  of  phases  of  several  prominent  point  scatterers  is  needed.  It 
can  be  accomplished  by  the  joint  time-frequency  analysis  [6].  With  the  detection 
algorithm,  we  have  the  ability  to  distinguish  the  time  intervals  when  the  target  undergoes 
smooth  2D  motion  from  those  containing  more  chaotic  3D  motion.  As  a  result,  the  good 
imaging  intervals  where  focused  images  are  more  easily  formed  can  be  automatically 
determined. 

The  paper  is  organized  as  follows.  First,  the  ISAR  imaging  problem  is  formulated  in 
terms  of  a  point  scatterer  model  in  Section  2.  In  Section  3,  the  2D  motion  assumption  in 
existing  motion  compensation  algorithms  is  analyzed.  We  show  the  reason  why  3D 
motion  is  a  problem  for  ISAR  imaging.  Section  4  discusses  the  3D  motion  detection 
algorithm  in  detail.  Examples  from  both  simulation  and  measurement  data  are  presented 
in  Section  5.  The  conclusions  are  given  in  the  last  section. 


2.  2D  and  3D  Motion  Models 

The  standard  model  used  in  ISAR  processing  is  the  point  scatterer  model  given  as 


E(f, to)  =  ^  ov(jr„^,)exp{  -j^—\r{tD)  +xt  +ynp(tD)\ j  (I) 

where /is  the  radar  frequency  and  to  is  the  dwell  time.  The  radar  echo  data  £(/  tD)  is 
in  the  (frequency /(dwell  time)  domain,  j c  and  y  represent  the  target  range  and  cross¬ 
range  positions,  respectively.  The  target  consists  of  Ns  point  scatterers,  with  the 
point  scatterer  depicted  by  position  (xh  yt)  and  strength  o{xh  yt).  The  target  motion 
includes  both  the  translational  motion  described  by  r(tD)  and  the  rotational  motion 
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Figure  I.  Geometry  of  an  JSAR  problem. 


described  by  When  there  is  no  translational  motion  and  the  rotational  motion  is 

uniform,  it  is  seen  that  a  2D  Fourier  transform  brings  the  radar  data  E(  f  tD)  into 
a  radar  image  o{x,  >>).  Otherwise,  motion  compensation  is  a  critical  step  in  ISAR 

imaging.  .... 

The  above  model  is  what  we  call  a  2D  problem  since  the  target  rotational  motion  is 
confined  to  a  2D  plane  and  describable  in  terms  of  only  one  angular  parameter  tp.  When 
there  is  3D  motion  of  the  target,  a  more  general  3D  model  is  required: 

E{f,  tD)  =  £  <xi(xi,yi,Zj)ex. p{  4- x-,  +  z/0(to)]}-  (2) 

/=  1 

In  the  above  expression,  a  third  coordinate  z  of  the  target  is  included  to  represent  the 
3D  target  and  another  independent  angular  motion  parameter  8  is  introduced  to 
describe  the  3D  rotational  motion  (see  Figure  1).  It  is  possible  to  perform  3D  target 
imaging  if  the  target  motion  is  known  exactly  [11],  [12].  In  practical  ISAR  scenarios, 
however,  we  have  no  access  to  the  target  motion.  Our  objectives  here  are  to  examine 
the  effect  of  3D  motion  on  ISAR  imaging  and  devise  an  algorithm  to  detect  the 
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Figure  2.  Illustration  of  2D  motion  vs.  3D  motion. 

it.  j  j  .Q  jtcptf  i  e  without  any  additional  knowl* 
presence  of  3D  motion  from  the  radar  data  itself,  t.e.,  wnnou  y 

edge  of  the  target  motion. 

3.  Problem  of  Existing  Motion  Compensation  Algorithms  with  3D  Target  Motion 

i.e., 

(3) 

0(tp)  =  by{tD)- 

This  allows  us  to  cast  equation  (2)  into  the  form 

the  3D  target  <r(x„  yh  z,)  onto  the  2D  motion  plane.  phase 

srr-r.LT-  * — -  -  *• » 

model  is  adequate.  _jel  appiicable  if  either  the 

rrr  srsr «=r=r  -  -*  -  *■ 


ALGORITHM  TO  DETECT  THE  PRESENCE  OF  3D  TARGET  MOTION 


227 


(a)  0» 


Figure  3.  Phase  linearizalion  achieved  by  rotational  motion  compensation,  (a)  Phases  of  two  point  scatterer  with 
2D  rotational  motion,  (b)  Both  phases  are  linearized  with  rotational  motion  compensation,  (c)  Phases  of  two  point 
scatterers  with  3D  motion,  (d)  Only  one  phase  is  linearized  with  rotational  motion  compensation. 


translational  motion  compensation  is  independent  of  the  models  in  (1)  and  (2),  only  the 
rotational  motion  compensation  needs  to  be  investigated.  With  2D  rotational  motion 
present,  the  phase  of  a  point  scatterer  i  due  to  the  rotational  motion  is 


Pi((d)  -yi'Pito)- 

Here,  the  constant  A-nf/c  has  been  suppressed  for  simplicity.  As  we  can  see  from  (5),  the 
phases  of  all  the  point  scatterers  are  linearly  related  to  each  other  (through  the  ratio  of 
their  cross  range  positions).  If  we  make  one  of  the  phases  a  linear  function  of  time,  then 
all  the  phases  are  linearized  simultaneously,  and  the  whole  target  can  be  focused  after  the 
Fourier  transform.  This  is  the  basis  of  most  2D  rotational  motion  compensation  algorithms 
based  on  the  point  scatterer  model  [3],  [4],  [5],  [6],  This  concept  is  illustrated  in  Figure  3. 
Figure  3(a)  shows  the  phase  functions  of  two  point  scatterers  under  2D  rotational  motion. 
Figure  3(b)  shows  that  both  points  can  be  made  linear  functions  of  time  after  we  force  one 
of  them  to  be  a  linear  function. 
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Figure  4.  Problem  with  a  typical  motion  compensation  algorithm,  (a)  Target  undergoes  a  2D  motion,  (b)  Image 
after  translational  motion  compensation,  (c)  Image  after  rotational  motion  compensation  (d)  Target  undergoes 
3D  motion,  (e)  Image  after  translational  motion  compensation,  (f)  Image  after  rotational  compensati  n. 


With  3D  motion,  the  phase  of  a  point  scatterer  due  to  the  rotational  motion  is 


Pi(*D)  =yiV(tD)  +Zi6{tD). 

In  this  case,  the  phases  of  the  point  scattered  are  no  longer  linearly  related.  If  we  make  one 
of  the  phases  a  linear  function  of  time,  the  phases  of  the  other  point  scattered  are  not 
automatically  made  linear  functions  of  time,  as  was  the  case  of  2D  motion.  Figure  3(c) 
shows  the  phase  functions  of  two  point  scatterers  with  3D  motion.  As  we  can  see  from 
Figure  3(d),  after  one  point  is  forced  to  be  of  linear  phase,  the  phase  of  the  other  point 

remains  nonlinear.  „  ,  .  ,  _ 

Figure  4  illustrates  some  simulation  results  of  the  effect  of  the  rotational  motion 
compensation  based  on  the  model  in  (1)  on  the  final  images  under  2D  mid  3D  target 
motion.  The  adaptive  joint  time-frequency  (AJTF)  algorithm  reported  in  [6]  ,s  used  for 
motion  compensation.  Ten  points  in  3D.  space  are  used  to  simulate  the  radar  data. 
Figure  4(a)  shows  an  assumed  2D  rotational  motion.  Figure  4(b)  shows  the  image  after 
the  translational  motion  compensation.  The  image  shows  one  point  being  focuse  in 
range  cell  25  while  other  points  are  unfocused  due  to  the  rotational  motion.  Figure  4(c) 
shows  the  image  after  the  2D  rotational  motion  compensation  in  which  a  point  scatterer 
in  range  cell  57  is  selected  for  focusing.  All  the  point  scatterers  are  focused  in  the 
image.  The  situation  with  an  assumed  3D  target  motion  is  shown  in  Figures  4(d)  (  ). 
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Figure  4(d)  shows  the  assumed  3D  motion.  Figure  4(e)  shows  the  image  after  trans¬ 
lational  motion  compensation.  Figure  4(f)  shows  the  final  image  after  the  2D  rotational 
motion  compensation.  The  two  points  in  range  cells  25  and  57  are  focused,  as  expected. 
Another  point  scatterer  in  range  cell  99  is  also  focused  as  it  happens  to  be  in  the  same 
2D  motion  plane  as  die  point  scatterer  in  range  cell  57.  However,  it  is  not  possible  to 
focus  all  the  points  simultaneously  with  an  existing  algorithm  based  on  the  2D  motion 
model. 


4.  3D  Motion  Detection  Algorithm 

Since  existing  motion  compensation  algorithms  cannot  handle  3D  target  motion,  it  is 
desirable  to  develop  a  general  compensation  algorithm  that  can  accommodate  3D 
motion.  However,  this  is  a  difficult  task  (see  [13]  for  background  on  this  problem) 
and  outside  the  scope  of  this  work.  Our  goal  here  is  to  develop  an  algorithm  to  detect  the 
presence  of  3D  motion  from  radar  data.  If  we  can  reliably  detect  those  time  intervals 
where  2D  target  motions  are  predominant,  we  can  use  the  existing  2D  motion 
compensation  algorithms  to  form  well-focused  ISAR  images. 

As  discussed  in  the  last  section,  2D  motion  can  be  represented  by  a  linear  relationship 
between  6  and  <£>.  Therefore,  we  set  out  to  detect  the  existence  of  a  nonlinear 
relationship  between  6  and  <p  in  our  3D  motion  detection  algorithm.  First,  we  write 
the  relationship  between  0  and  ip  into  a  linear  and  a  nonlinear  part  as  follows: 

0(to)  =  bipito)  +  m(to)  (?) 

where  b  is  the  linear  constant  and  m(tD)  is  the  nonlinear  part  which  indicates  deviation 
from  2D  target  motion,  or  the  degree  of  3D  motion.  Next  we  try  to  gather  target  motion 
information  by  analyzing  the  phases  of  two  point  scatterers  on  the  target.  Let  us  write 
the  relationship  between  the  phase  functions  P\  and  P2  of  two  point  scatterers  as: 


Pi(to)  =  aP\  (*d)  +  (*0 

The  relationship  is  again  decomposed  into  the  linear  part,  where  a  is  the  linear  constant, 
and  the  nonlinear  part  n(tD).  Our  goal  is  to  derive  a  relationship  between  m(tD)  and  n(tD) 
so  that  the  presence  of  m  can  be  detected  by  observing  n. 

After  the  standard  translational  motion  compensation,  the  time-varying  phase  of  a  point 
scatter  is  in  the  form  of 


Pi(tD)  =  A ynp{tD)  +  A  Zi&(tD)  (9) 

where  Ay,  and  Az,  are  differential  positions  of  point  scatterer  i  relative  to  the  reference 
point  chosen  during  translational  motion  compensation.  Substituting  (7)  into  (9)  and  then 
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evaluating  (9)  at  point  scatterers  1  and  2,  we  have 

P\(tD)  =  (Aji  +  bAzi)<p(tD)  +  Az\m{tD)  (10a) 

P2{tD)  =  (Ay2  +  bAz2)p{tD)  +  A  z2m(tD).  ( ,ob) 

We  next  substitute  (10)  into  (8),  which  leads  to 

a[(Ay\  +  bAz\)y{tD)  +  A  z,m(f0)]  +  «(*o)  =  (Ay2  +  Mz2)y?(to)  +  Az2  m{tD)  (1 1) 

Notice  that  if  there  is  only  2D  motion,  then  the  phases  of  the  two  point  scatterers  must  be 
linear.  This  means  if  m  =  0,  then  n  =  0.  By  using  this  fact  and  equating  the  coefficients  of 
ip(tD)  in  (11),  the  constant  a  can  be  derived: 

_  Ay2  +  bAz2  (12) 

°  Ay\  +  bAz\  ' 

By  substituting  (12)  into  (11),  we  finally  arrive  at 

Ayi+frAzi  .  (13) 

m  Az2Ay\  -  AyjAz\ 

Equation  (13)  states  that  once  the  nonlinear  phase  term  n  is  known,  it  is  proportional  to 
nonlinear  motion  m.  Therefore,  the  steps  to  determine  the  degree  of  3D  target  motion  are 
as  follows.  First,  we  extract  the  phases  of  two  point  scatterers  from  the  radar  data.  Next  we 
find  the  nonlinear  phase  function  n  using  a  minimum  least  squares  fit  of  equation  (8).  Once 
n  is  known,  we  use  equation  (13)  to  decide  on  the  degree  of  3D  motion.  The  remaining 
issues  are:  (i)  how  to  determine  the  phase  functions  of  the  point  scatterers,  (ii)  how  to 
define  the  degree  of  nonlinearity  and  the  degree  of  3D  motion  once  n  is  known,  and  (iii) 
how  to  compare  the  degree  of  3D  motion  from  one  imaging  interval  to  another.  These 
three  issues  are  discussed  in  the  following  subsections. 

4.1.  Phase  Estimation  Using  Adaptive  Joint  Time-Frequency  Projection 

After  the  translational  motion  compensation,  the  radar  signal  contains  only  rotational 
motion.  To  estimate  the  phase  of  a  prominent  point  scatterer,  we  utilize  the  adaptive  joint 
time-frequency  (AJTF)  projection  technique  discussed  in  [6].  We  begin  with  the  radar  data 
in  the  (range)-(dwell  time)  domain.  Within  a  fixed  range  cell,  the  data  can  be  written  as 

£,(to)  =  &i  exp  (  ~J  ~~  (yiVito)  +  z/0(4>))  ( 1 

1=1 

where  fc  is  the  center  frequency.  Among  the  Nr  point  scatterers  within  the  range  cell. 
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Figure  5.  (a)  (Dwell  trmeMDoppler  frequency  representation  of  radar  signal  in  a  range  cell  with  three  point 
scatterers.  (b)  The  basis  function  that  is  best  matched  to  the  dominant  point  scatterer  is  found  by  the  AJTF  project 
method. 

we  express  the  phase  behavior  of  the  strongest  one  as  a  polynomial  function: 

</>mUd)  =  (fl*  +fltZ  +fi^  +  ••) 
and  consider 

h{t)  =  exp(-_/V>A,  (/£>)] 

as  a  basis  for  the  radar  signal.  The  phase  parameters  are  then  found  by  searching  for  the 
maximum  projection  from  the  radar  signal  onto  the  basis  function. 


ifiifii . . .)  =  arg  max  J  E{tD)h* {tD)dtD  . 


Equation  (16)  means  that  the  phase  iunction  parameters  are  estimated  to  give  a  maximum 
projection  from  the  radar  data  onto  the  basis  function  for  that  prominent  point  scatterer.  In 
die  search  procedure,  the  first  term  /,  can  be  obtained  by  using  the  fast  Fourier  transform, 
while  all  other  higher  order  terms/2,/3, ...  are  obtained  using  exhaustive  search.  Figure  5 
illustrates  the  process  of  AJTF  phase  estimation.  Figure  5(a)  shows  the  radar  signal  in  one 
range  cell  with  three  point  scatterers  in  the  joint  (dwell  time)-(Doppler  frequency)  plane. 
The  tilted  trajectory  of  the  prominent  point  scatterer  1  implies  there  exist  higher-order 
terms  in  the  phase  function.  Figure  5(b)  shows  the  basis  fiinction  h(tD).  During  the  search, 
we  change  the  position  (/,),  tilting  (/2)  and  curvature  (/3, . . .  )  of  h  until  the  projection  of 
h  onto  the  radar  signal  is  maximized. 

4.2.  Measure  of  Nonlinearity  Between  Two  Phase  Functions 


We  notice  that  in  (8),  the  two  phase  functions  are  formulated  with  a  linear  relationship  plus  a 
nonlinear  residual  part.  After  the  two  phase  functions  are  estimated  using  the  AJTF  tech¬ 
nique,  a  least-squares  fitting  can  be  performed  to  generate  die  best-fit  linear  part.  The  actual 
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Figure  6.  Measure  of  the  nonlinearity  of  two  phase  functions. 

phases  deviates  from  this  linear  relationship.  The  deviation  n  is  inte^ted  over  tte  dwell 
time  to  represent  the  degree  of  phase  nonlinearity  over  the  imagmg  interval  as  follows. 

N\i  —  J  \nn(tD)\dtD-  ^  ^ 

The  process  is  illustrated  in  Figure  6.  The  solid  line  is  the  actual  relationship 

phase  functions  P\  and  P2.  The  dotted  line  is  the  linear  approximation  of  the  relationship. 

The  area  of  the  shadowed  region  is  A/)  2* 

In  a  similar  fashion,  we  define  the  degree  of  3D  motion  as  the  deviation  from  a  linear 
relationship  between  6  and  over  the  dwell  interval  as  follows: 

M=  J  \m(tD)\dtD ■  ^ 

Based  on  (13),  we  see  that  M  and  re  directly  related. 

Nn  =  PnM  (19a) 


where 

_  AziAy2  -  Avi^22  (19b) 

^I2  ~  Ayi  +  bkz\ 

Thus  by  finding  the  observable  Ni2,  we  can  obtain  the  degree  of  3D  motion  M  to  within  a 
proportionality  constant. 
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Figure  7.  Effect  of  the  number  of  point  scatterers  used  on  3D  motion  detection  result. 

4.3.  3D  Motion  Comparison  among  Different  Imaging  Intervals 

As  indicated  by  (19),  the  phase  nonlinearity  of  two  point  scatterers  N  is  proportional  to  the 
degree  of  3D  motion  M,  so  we  can  use  the  detected  phase  nonlinearity  as  a  measure  of  3D 
motion.  However,  we  notice  that  the  constant  of  proportionality  is  dependent  on  the  point 
scatterer  positions.  A  problem  arises  when  we  need  to  compare  the  detection  result  from 
one  imaging  interval  to  that  from  another  imaging  interval.  Since  we  cannot  guarantee  that 
we  track  the  same  set  of  points  from  frame  to  frame,  the  proportionality  constant  can 
change  from  frame  to  frame,  and  we  cannot  reliably  observe  M  from  W  across  femes.  To 
overcome  this  difficulty,  we  track  more  than  two  point  scatterers  within  each  feme  and 
compute  Ny  for  each  pairing  of  scatterers  i  and  j  (i  ^  j).  Then  we  generate  an  average 
value  <Ny>  from  all  the  possible  phase  relationships.  From  (19),  we  have 


(Ny)  =  (0ij)M. 

We  postulate  that,  from  a  statistical  point  of  view,  </%>  approaches  a  constant  that  is 
independent  of  frames  if  we  average  over  a  sufficient  number  of  point  scatterers.  If  this  is 

true,  <Ny>  should  become  a  good  indicator  of  M. 

We  test  the  effectiveness  of  this  approach  on  the  detection  result  by  simulation.  We  input  a 
set  of  motion  parameters  and  generate  the  phase  functions  based  on  the  3D  motion  model 
20  point  scatterers  from  an  airplane  model  is  used.  We  then  randomly  choose  a  number  of 
point  scatterers  and  use  their  phase  functions  to  compute  <NU>.  We  examine  how  the 
results  vary  as  different  number  of  point  scatterers  is  used.  We  find  that  the  results  begin  to 
converge  after  about  5  scatterers.  Figure  7  shows  a  plot  of  <N0>  versus  the  frame  number  if 
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Figure  8.  Data  processing  flow  chart. 


we  use  5  point  scattered  ( 1 0  phase  pairs).  If  we  increase  the  number  ofpomt  scatterersto  10 
(45  phase  pairs),  there  is  only  minor  change  in  the  detection  output.  Therefore,  <N0> 
he  used  to  indicate  the  degree  of  3D  motion  given  a  sufficient  number  of  point  scatterers. 


(a) 


Figure  9.  (a)  Detected  3D  motion  from  simulated  radar  data,  (b)  Degree  of  3D  motion  from  truth  motion  data. 
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5.  Results 

To  demonstrate  the  effectiveness  of  the  3D  motion  detection  algorithm,  we  test  our  algorithm 
on  radar  data  from  two  targets.  The  first  target  is  an  aircraft,  which  flew  in  a  large  clockwise 
circle  during  a  9-minute  interval.  We  also  have  access  to  the  target  motion  data  through  the 
GPS  (global  positioning  system)  and  INS  (inertial  navigation  system)  sensors  carried  on¬ 
board  the  aircraft  [14].  Figure  8  shows  our  processing  flow  chart.  The  GPS/INS  data  is  used 
to  establish  the  truth  target  motion.  The  raw  radar  data  is  used  as  input  to  the  3D  motion 
detection  algorithm.  We  can  also  generate  the  1SAR  images  using  our  AJTF  motion 
compensation  algorithm  [6].  We  are  therefore  able  to  both  compare  the  detection  result 
with  the  truth  motion,  and  observe  the  effect  of  the  3D  motion  on  the  ISAR  image  quality. 

We  first  test  the  3D  motion  detection  algorithm  on  simulated  radar  data.  To  generate  the 
simulation  data,  we  use  the  actual  motion  data  from  the  GPS/INS  sensors  in  conjunction 
with  a  point  scatterer  model.  From  the  aircraft  model,  60  point  scatterers  are  selected  to 
simulate  the  radar  data  based  on  the  actual  motion  data  and  equation  (2).  Five  range  cells  are 
then  chosen  for  phase  analysis  in  the  detection  procedure.  Figure  9(a)  shows  the  detected 
degree  of  3D  motion  for  20  image  frames  from  the  simulated  radar  data.  For  comparison. 
Figure  9(b)  shows  the  degree  of  3D  motion  obtained  based  on  the  truth  motion  date.  The 
frames  with  significant  3D  motion  are  highlighted  with  circles  and  the  frames  with  2D 
motion  are  highlighted  with  diamonds.  It  is  seen  that  the  two  results  agree  fairly  well. 

Next,  we  test  the  detection  algorithm  using  the  actual  radar  measurement  data. 
Figure  10(a)  shows  the  detected  3D  motion  from  the  radar  data  over  20  frames.  The 
corresponding  imaging  interval  for  each  frame  is  2.3  seconds  while  the  total  flight  duration 
is  5  minutes.  The  four  frames  with  the  most  significant  3D  motion  based  on  our  detection 
algorithm  are  labeled  as  circles.  They  are  frames  6,  14, 17  and  18.  Figure  10(b)  shows  the 
degree  of  3D  motion  obtained  based  on  the  truth  motion  data.  We  observe  that  the  truth 
motion  date  indeed  contains  a  high  degree  of  3D  motion  at  those  four  frames  detected  by 
our  algorithm. 

To  further  examine  the  quality  of  the  ISAR  images  when  3D  motion  is  present,  we 
generate  images  using  our  motion  compensation  algorithm  in  Figures  11  to  14.  Figure 
1 1(a)  shows  the  plot  of  6  vs.  <f)  derived  from  the  truth  motion  date  for  frame  1 8,  which  is  a 
frame  found  to  contain  substantial  3D  motion.  The  actual  motion  is  shown  in  the  solid 
curve  and  the  dashed  line  is  the  best-fit  2D  motion  approximation.  It  is  clear  that  the  curve 
deviates  significantly  from  the  dashed  line  and  the  actual  motion  cannot  be  well 
approximated  with  2D  motion.  Figure  11(b)  shows  the  resulting  image  obtained  after 
the  motion  compensation,  and  is  blurred  in  the  Doppler  dimension  (vertical  axis).  As 
expected,  die  2D  motion  compensation  algorithm  cannot  focus  all  the  points  due  to  the  3D 
target  motion.  Figures  12(a)  and  12(b)  show  the  same  conclusion  for  frame  14,  which  is 
another  frame  identified  as  having  significant  3D  motion.  In  Figure  13,  we  show  the 
results  for  frame  2,  which  has  very  little  3D  motion.  As  we  can  see  from  Figure  13(a),  the 
actual  motion  can  be  well  approximated  by  a  line  in  the  Q-p  plot.  The  image  shown  Figure 
13(b)  is  well  focused.  In  particular,  the  point  scatterers  on  the  target  show  nearly  equal 
range  and  Doppler  extent,  contrary  to  the  previous  two  images.  The  aircraft  body  line  is 
clearly  recognizable. 
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Figure  12.  3D  motion  and  resulting  ISAR  image  (frame  14). 


From  Figure  10,  we  notice  that  there  exists  a  discrepancy  in  frame  11,  where  the 
detection  result  does  not  indicate  any  3D  motion  while  the  truth  motion  data  shows  a 
significant  amount  of  3D  motion.  The  truth  motion  is  shown  in  Figure  14(a),  confirming 
the  presence  of  3D  motion.  One  explanation  is  that  those  prominent  points  used  by  the 
detection  algorithm  lie  nearly  on  a  2D  plane  so  that  they  still  can  be  focused.  As  we  have 
discussed  in  Section  3,  the  2D  model  is  applicable  if  either  the  motion  is  2D  or  the  target  is 
of  2D  in  extent.  It  is  likely  that  the  latter  condition  is  met  for  this  frame.  This  is  confirmed 
by  the  image  shown  in  Figure  14(b).  We  see  that  the  image  quality  is  actually  not  so  bod. 
Therefore,  our  detection  algorithm  objectively  reflects  the  quality  of  the  images  generated 

bv  the  2D  motion  compensation.  .  . 

A  second  data  set  is  used  to  test  our  3D  motion  detection  algorithm.  This  data  set 
consists  of  the  ISAR  data  collected  from  a  small  ship  on  the  ocean.  Because  of  the  surface 
movement  of  the  sea,  the  target  is  believed  to  have  considerable  3D  motion  during  the 
imaging  intervals.  The  3D  motion  detection  result  is  shown  in  Figure  15(a)  with  the  peaks 


(a)  (b) 


Figure  13.  2D  motion  and  resulting  ISAR  image  (frame  2). 


239 


ALGORITHM  TO  DETECT  THE  PRESENCE  OF  3D  TARGET  MOTION 

con-esponding  to  regions  with  3D  motion.  The  total  data  duration  is  20  seconds  and  the 
imaging  dwell  time  is  0.64  second  per  frame.  For  this  data  set,  reliable  truth  target  motion 
is  not  available.  Instead,  we  generate  the  motion  compensated  images  shown  m .Figures 
15(b)  to  15(d)  to  demonstrate  the  effect  of  3D  motion  on  ISAR  image  quality.  The  image 
frame  with  the  largest  detected  3D  motion,  frames  3,  is  shown  m  Figure  15(b).  It  is  poor  y 
focused  Figure  15(c)  shows  the  image  from  frame  14,  which  is  the  frame  with  the  secon 
highest  detected  3D  motion.  The  frame  with  the  smallest  3D  motion  based  on  our 
algorithm,  frame  20,  is  shown  in  Figure  15(d).  It  shows  a  well-focused  ISAR  image.  This 
left  confirms  the  effectiveness  of  our  algorithm  in  detecting  good  imaging  intervals  from 
those  imaging  intervals  containing  large  3D  motion. 


6 .  Conclusions 

In  this  paper,  we  set  out  to  develop  an  algorithm  to  detect  the  presence  of  3D  target  motion 
from  ISAR  data.  Based  on  the  3D  point  scatterer  model,  we  first  examined  the  effect  of  3D 
motion  on  ISAR  imaging.  It  was  shown  that  the  existing  motion  compensation  algorithms 
could  not  properly  focus  targets  exhibiting  3D  motion  during  the  imaging  interval.  We 
then  derived  an  algorithm  to  blindly  detect  the  degree  of  3D  target  motion  from  raw  radar 
data.  It  is  based  on  measuring  the  linearity  of  phases  between  two  or  more  point  scatterers 
on  the  target.  The  phase  estimation  was  implemented  using  the  adaptive  joint  time- 
frequency  technique.  Examples  were  provided  to  demonstrate  the  effectiveness  of  the  3 
morion  detection  algorithm  with  both  simulation  and  real  ISAR  data.  The  detection  results 
were  corroborated  with  the  truth  motion  data  from  on-board  motion  sensors  and  coirclated 
with  the  resulting  ISAR  images.  With  the  detection  algorithm,  we  have  the  ability  to 
distinguish  the  time  intervals  when  the  target  undergoes  smooth  2D  motion  from  those 
containing  more  chaotic  3D  motion.  As  a  result,  the  good  imaging  intervals  where  focused 
images  are  more  easily  formed  can  be  automatically  selected. 
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Use  of  Genetic  Algorithms  in  ISAR  Imaging  of 
Targets  with  Higher  Order  Motions 


t7rnrt»T  dgoriflro  (GA)  me  proposed  for  inverse  synthetic 
aperture  radar  (ISAR)  imaging.  GA  h  osed  for  matin 
parameters  search  in  place  rf  eahauath*  search  in  the  adaptive 
jaint  ttaMTeqaeney  (AJTF)  algorithm.  While  ajaint^iug  the 
same  accuracy,  G  A  has  lower  computational  complexity,  especially 
for  targets  with  higher  ssder  motions. 


I.  INTRODUCTION 

An  Inverse  synthetic  aperture  radar  (ISAR)  system 

usually  collects  radar  data  from  a  target  moving  on 
the  ground,  in  the  air,  or  over  the  ocean.  In  the  ISAR 
problem,  the  radar  is  stationary  while  the  target  moves 
with  both  translation  motion  and  rotational  motion. 

In  the  microwave  frequency  range,  ISAR  has  been 
identified  as  an  effective  tool  for  target  identification 

HI 

Target  motion  is  an  essential  part  in  ISAR 
imaging.  On  die  one  hand,  target  motion  is  needed 
to  generate  Doppler  (or  cross-range)  resolution  [2]. 

On  the  other  hand,  unwanted  motion  causes  image 
blurring.  When  the  target  has  uniform  rotational 
motion  only  and  the  radar  data  is  collected  over  a 
small  angular  aperture,  a  simple  Fourier  transform 
will  bring  the  raw  radar  data  into  a  two-dimensional 
ISAR  image.  However,  actual  targets  observed  by 
operational  radar  rarely  have  such  an  ideal  motion. 
Therefore,  motion  compensation  is  needed  to  generate 
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focused  ISAR  imagery.  There  exist  many  different 
motion  compensation  algorithms  to  deal  with  target 
motion  [3-5].  Most  of  the  algorithms  start  with  coarse 
range  alignment  based  on  the  correlation  of  the  range 
profiles.  Then  die  phase  information  within  one  range 
cell  is  utilized  to  achieve  fine,  motion  compensation. 

Phase  estimation  is-  critical  in  fine  motion 
compensation.  Compared  with  die  amplitude;  the  . 
phase  of  the  radar  signal  is  much  more  sensitive  to 
die  change  in  range.  Based  on  the  concept  of  signal 
parameterization  reported  in  [6, 7],  an  adaptive  joint 
time-frequency  (AJTF)  .algorithm  was  proposed 
in  [8]  for  phase  estimation  of  a  prominent  point 
scarterer.  In  this  method,  the  target  motion  is  modeled 
as  a  polynomial  function  and  an  exhaustive  search 
procedure  is  used  to  find  the  motion  parameters  that 
are  embedded  in  the  phase  of  the  prominent  point 
scatterer.  While  the  performance  of  this  algorithm  is 
very  good  [9],  the  main  bottleneck  in  this  procedure 
is  the  computational  complexity  associated  with  the 
parameter  search.  When  die  target  motion  is- highly 
irregular,  i.e.,  when  the  number  of  parameters  needed 
to  model  the  motion  is  large,  the  use  of  die  exhaustive 
search  becomes  prohibitively  expensive. 

Our  objective  here  is  to  reduce  die  computation 
time  associated  with  the  motion  parameter  search  , 
in  the  AJTF  procedure.  Our  proposed  approach  is 
to  incorporate  genetic  algorithms  (GA)  [10]  mto 
the  AJTF  search  process.  (Some  preliminary  work 
on  this  topic  was  reported  in  [1 1].)  In  contrast 
with  conventional  -optimization  methods,  GA  is  a 
population-based,  statistical  search  technique.  It  . 
borrows  such  concepts  as  inheritance  and  mutation 
from  the  biological  .evolution  process  .[12],  As  a 
global  optimization  technique,  GA  is  known  to  be  ; 
very  easy  to  implement  and  applicable  to  many  design 

and  inverse  problems  [13],  . 

This  paper  is  organized  as  follows.  In  Sections 
n  and  III,  we  outline  the  methodology.  The  AJTF 
analysis  for  ISAR  motion  compensation  is  described 
in  Section  1L  GAs  are  introduced  in  Section  m. 

The  next  two  sections  include  results  and  analysis. 

In  Section  IV,  simulations' with  point  scatterers 
are  provided  to  validate  the  use  of  GA  for  phase 
estimation  Measurement  data  processing  results  are 
shown  in  Section  V.  Finally,  conclusions  are  given  in 
Section  VI. 

II.  ISAR  MOTION  COMPENSATION  USING  jOlNT 
TIME-FREQUENCY  PROJECTION 

We  restrict  our  attention  here  to  rigid  body  targets 
with  a  fixed  rotational  axis.  We  also  assume  that 
the  target  undergoes  only  a  small,  angular  rotation 
within  the  dwell  interval.  Under  these  assumptions. 


lWe  consider  the  conventional  motion  assumptions  for  ISAR 
imaging  in  this  paper.  Those  assumptions  may  not  hold  ^  rawing 
Homtd  vehicles  or  small  ships.  The  reader  is  referred  to  [18-20)  for 
a  detailed  discussion  of  those  more  challenging  scenarios. 


a  two-dimensional  point  scatterer  model  relates  the 
radar  data  to  a  moving  target  through  the  following 
equation:  * 

E(f,tD)  -  XX  exp  |  -j~-[r(tD)  + 

i  *• 

(D 

where  /  is  the  frequency  and  tD  is  the  dwell  time.  In 
this  model,  the  radar  data  is  comprised  of  the  sum  of 
responses  from  a  collection  of  point  scatterers.  (Xj.yt) 
denotes  the  point  scatterer  position  while  <r,  denotes 
the  scatterer  strength.  The  target  motion  includes  both 
translation  motion  r(tD)  and  rotational  motion  y<rD). 

After  range  compression  and  range  alignment  to 
place  all  the  point  scatterers  in  their  correct  range 
bins,  the  radar  signal  through  one  range  cell  r  can  be 
expressed  in  die  form  of 

Er(tD)  =  £X;exp{— j^k[Ar(tD)  +  xt  +  y/V^tp)]  j 

(2) 

where  f0  is  the  center  frequency  of  the  radar  and 
the  index  includes  only  those  point  scatterers  in  the 
particular  range  celL  The  residual  translation  motion 
is  depicted  as  Ar(rD).  After  such  coarse  alignment 
procedure,  tbe  residual  translation  motion  is  smaller 
than  the  range  resolution.  However,  it  can  still  be 
larger  than  the  radar  wavelength.  Both  the  residual 
translation  motion  A r(tD)  and  the  rotational  motion' 
tp(f  D)  can  be  expanded  into  polynomial  functions  of 
the  dwell  time  as 

A r(tD)  =  attD  +  Ot?d  +  a3ip ...  ^ 

=  b\lt>  +  Vo  +  h fD  •  •  • 

where  any  coefficients  beyond  die  first  order  are 
detrimental  to  ISAR  image  formation.  To  solve  the 
ISAR  motion  compensation  problem,  we  need  to 
determine  these  motion  parameters  and  to  remove  the 
unwanted  nonlinear  phase  terms  from  the  radar  data. 

This  task  can  be  accomplished  using  the  AJTF 
procedure  [8]-  The  essential  idea  of  this  procedure 
is  to  find  the  basis  function  that  most  resembles  the 
strongest  signal  component  in  (2).  For  our  problem,  a 
basis  function  in  the*  form  of 

h(tD)  =  exp  |-j— ~ +  *  •  -)j  W 

is  used.  The  best  basis  is  found  by  searching  for 
parameters  /,,  /2,  /3v~that  maximize  the  projection 
from  the  radar  signal  onto  the  basis,  i.e., 

(/i»/2,/3,  -•>  =  aigmaxKEf(^),hM|  (5) 

where  the  projection  is  formulated  as  the  inner 
product  of  the  two  functions  as 

(Er(tD)MtD))  =  J  Er(tDW(tD)dtD.  •  (6) 
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Fig.  J.  Flowchart  of  GA 


In  (5),  the  linear  coefficient  /,  can  be  found  efficiently 
with  die  fast  Fourier  transform  (FFT).  Coefficients 
for  the  nonlinear  phase  terms,  f2,  .&>•••»  must  be 
determined  through  a  more  time-consuming  search. 

After  the  motion  parameters  of  the  prominent  point 
scatterer  are  found,  we  can  carry  out  the  translation 
motion  compensation  by  multiplying  the  radar  data 
with  the  conjugate  of  this  basis.  Since  all  the  point 
scattered  share  the  same  translation  motion  in  (2),  this 
operation  will  remove  the  translation  motion  of  the 
whole  target  Rotational  motion  compensation  can  also 
be  carried  out  by  estimating  the  phase  of  a  second 
point  scatterer  in  a  different  range  cell.  After  the  phase 
estimation,  we  can  resample  the  data  in  dwell  time  to 
make  the  phase  linear  [3, 4]. 

tit.  USE  OF  GENETIC  ALGORITHMS  FOR  PHASE 

PARAMETER  SEARCH 


ItiJtlOl  -»  fOfllltl  Mrtttion  c  -t  «P  +  (!-o)P. 

(a)  (b) 

R®.  2.  Examples  of  crossover  and  imitation  operations  in  binary 
and  real-coded  GA,  (a)  Binary-coded  GA.  (b)  Real-coded  GA. 

such  as  entropy  can  also  be  used).  The  search 
is  carried  out  in  a  range  cell  with  a  prominent 
point  scatterer.  When  such  a  good  range  cell  is 
not  available,  multiple  range  cells  can  be  used  to 
improve  the  GA  phase  estimation  accuracy.  If  we 
are  satisfied  with  the  solution,  the  process  is  done. 
Otherwise,  a  new  generation  is  produced  for  the 
next  evaluation-reproduction  iteration.  To  make 
op  individuals  in  the  new  generation,  some  good 
parents  are  selected  from  the  previous  generation 
and  two  operations  called  crossover  and  mutation  are 
applied  to  produce  children  in  the  next  generation. 
Whether  or  not  crossover  and  mutation  occur  is 
determined  randomly.  The  crossover  and  mutation 
probabilities  are  chosen  based  on  the  tradeoff  between 
two  conflicting  requirements.  Increasing  the  variation 
in  the  new  generation  brings  a  chance  for  better 
solutions,  but  it  tends  to  lose  the  features  of  the  good 
solutions  from  the  .previous  generation. 

Roughly  speaking’  there  are  two  kinds  of  GA. 

One  is  binary-coded  GA  [10}.  The  other  is  real-coded 
GA  [15].  In  the  former,  the  physical  parameters  to  be 
searched  are  first  discretized  into  binary  bits.  There  is 
a  one-to-one  mapping  between  a  physical  parameter 


As  we  have  discussed  in  the  last  section; 

ISAR  motion  compensation  can  be  formulated 
as  a  parameterization  process  for  both  translation 
motion  and  rotational  motion.  In  [8],  a  brute-force  : 
search  procedure  is  employed  to  carry  out  the  * 
parameterization.  This  means  that-  we  exhaustively  ^ 
search  the  solution  space  for  the  maximum  projection. 
G As* are  investigated  here  to  search  for  the  motion  .  - 
parameters  to  reduce  the  computation  time.  We  should 
point  out  that  although  a  structured*  tree-search  is 
an  easy  and  straightforward  way  to  decrease  the 
computational  complexity,  it  does  not  always  lead'to 
a  global  optimum.  ' 

GA  is  a  global  optimization  method  based  on  - 
concepts  from  ecological  systems  [10,  12, 14].  The 
flowchart  of  a  typical  GA  process  is.  Illustrated  in 
Fig.  1.  It  starts  by  setting  up  the  parameters  for  both 
the  physical  problem  ami  die  GA  implementation. 

GA  operates  on  a  population  of  many  individuals. 

The  initial  population  is  randomly  generated  within 
the  searching  space.  The  goodness  of  the  solution 
is  then  evaluated  based  on  an  objective  function.  < 
The  objective  function  is  chosen  as  the  projection  *  * 
magnitude  from  the  radar  signal  onto  the  basis 
function  (see  (5»  (althou^i  other  objective  functions 


C  and  its  AT-bit  binary  representation  c,,c2, as  . 
follows:  . 


i=i 


where  die  is  the  valid  search  space  for  C. 

A  candidate  solution  of  the  problem  is  expressed  in 
the  form  of  a  chromosome, k  the  collection  of 
bits  representing'all  the  parameters.  For  crosk>ver,  a 
random  break  point  in  the  chromosome  is  picked.  The 
bits  before  the  point  from  one  parent  are  combined  ^ 
with  the  bits  after  the  point  from  another  parent  to 
form  one  child.  Another  child  is  generated  in  the 
reverse  fashion.  For  mutation,  a  single  bit  is  randomly 
picked  and  its  valued  is  inverted.  The  crossover  and 
mutation  operations  for  binary-coded  GA  are  depicted 

in  Fig.  2(a).  ..  .  .  , 

In  real-coded  GA,  there  is  no  coding-decoding 
process  and  the  algorithm  directly  operates  on 
the  physical  parameters.  For  crossover,  a  linear 

combination  is'usiially  used  as  follows 

*  /  * 


•  ■  •Cl'^aP1  +{l-a)f^ 
C^U-ajPj  +  aP, 


(8) 
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fig.  3.  GA  for  third-order  ph«-  estimation,  (a)  Multi-modal  objective  function,  (b)  GA  convergence  curves,  (c)  GA-estimated  phase 
B  compared  to  the  original  truth  phase. 


where  two  children  (C,,C2)  are  reproduced  from  two 
parents  (PVP2).  a  is  a  random  number  between  0  aiid 
1  to  ensure  that  the  new  parameters  will  not  fall  out 
of  range.  For  mutation,  a  child  C  that  is  different  from 
the  parent  P  is  needed.  For  this  purpose,  a  solution  P, 
is  picked  up  randomly  from  the  searching  space  and 
linearly  combined  with  P  to  generate  C  in  the  same 
way  as  described  by  (8).  The  crossover  and  mutation 
operations  for  real-coded  GA  are  depicted  in  Fig.  . 
2(b). 

The  basic  theory  in  GA,  the  schemata  theory, 
seems  to  favor  the  use  of  binary-coded  GA 
[14].  Most  work  on  GA  has  followed  this  path. 
Recently,  researchers  have  also  experimented  with 
real-coded  GA  and  have  observed  some  advantages 
in  convenience  and  accuracy  [16],  In  the  next  section, 
we  test  both  binary-coded  and  real-coded  GA  in  our 
phase  estimation  problem. 

The  GA  process  is  usually  stopped  using  criteria 
based  on  the  performance  of  the  available  solutions  in 
the  present  generation.  In  our  case,  we  do  not  know 
the  true  maximum  projection  value.  Therefore,  we 
choose  to  stop  the  GA  process  when  the  projection 
value  does  not  increase  after  a  certain  number  of 
generations. 


IV.  POINT  SCATTERER  SIMULATION 

Point  scatterer  simulations  are  first  used  to  test 
the  use  of  GA  for  ISAR  motion  compensation.  We 
first  test  the  accuracy  of  GA  phase  estimation.  In  this 
example,  we  use  two  point  scatterers  with  amplitudes 
1  and  02.  They  are  located  within  one  range  cell 
and  contain  third-order  translation  motion  (i.e.,  the 
coefficients  a,,  and  03  in  (3)  are  significant  while 
all  higher  order  coefficients  are  zero).  We  run  both 
binary-coded  GA  and  real-coded  GA  to  search  for 
Qj  and  for  this  simple  phase  estimation  problem. 
The  population  size  is  50.  In  both  cases,  the  crossover 
probability  is  0.8  and  the  imitation  probability  is  0.3. 
For  exhaustive  search,  the  search  for  and  a3  is 
carried  out  on  a  128  by  128  grid.  We  choose  a  7-bit 
representation  in  binary-coded  GA  and  search  in  the 
same  discrete  space  as  in  the  exhaustive  search.  The 
actual  objective  function  surface  is  plotted  in  Fig.  3(a). 
We  observe  many  local  maxima,  indicating  this  would 
be  a  challenging  problem  for  a  local  optimization 
method.  Fig.  3(b)  shows  the  GA  convergence  curve, 
with  the  real-coded  GA  result  in  solid  line  and  the 
binary-coded  GA  result  in  dashed  line.  Since  GA  is 
a  statistical  approach,  we  do  not  get  exactly  the  same 
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Fig.  4.  Performance  of  GA  compared  to  exhaustive  search,  (a)  Computational  complexity  as  a  function  of  the  number  of  parameters. 

(b)  Accuracy. 


result  from  each  GA  run.  To  decrease  the  statistical 
variation,  the  results  in  Fig.  3(b)  are  obtained  by 
averaging  over  20  GA  runs.  We  can  see  that  after 
about  150  generations  the  two  projection  corves 
nearly  converge  to  the  truth  value  of  1.  We  also 
observe  that  real-coded  GA  produces  a  slightly  higher 
projection  value.  Fig.  3(c)  shows  the  resulting  phase 
from  a  single  GA  ran  after  200  generations.  The 
estimated  phase  from  binary-coded  GA  is  plotted 
in  circles,  the  phase  from  real-coded  G A  in  crosses, 
and  the  original  phase  function  in  solid  line.  We  see 
very  good  agreement  among  the  three  results,  meaning 
good  accuracy  from  the  two  GA  results. 

,  In  the  second  example  we  compare  the 
computational  complexity  of  GA  to  exhaustive  search 
for  different  orders  of  motion.  As  we  have  pointed 
earlier,  the  main  problem  with  exhaustive  search 
for  motion  parameter  extraction  is  the  computational 
load.  This  problem  becomes  more  acute  when  the 
order  of  the  motion  is  high.  Again,  we  use  two  point 
scatterers  with  amplitudes  1  and  0.2.  We  generate 
the  radar  data  from  these  two  point  scatterers  with 
some  preset  motion.  We  then  apply  exhaustive 
search,  binary-coded  GA,  and  real-coded  GA  for  the 
phase  estimation  problem  with  different  orders  of 
motion.  The  same  GA  parameters  as  in  the  previous 
example  are  used  and  the  results  are  averaged  over 
20  runs.  The  exhaustive  search  is  known  to  have  an 
exponential  complexity  of  f?(exp(n)).  As  expected, 
the  resulting  computation  time  in  logarithm  scale 
shews  up  as  a  straight  line  in  Fig.  4(a).  For  GA,  no 
theoretical  complexity  formulation  is  available  in 
general.  (The  complexity  of  0(«logn)  is  claimed 
for  selected  test  functions  in  [17].)  We  run  both  - 
binary-coded  GA  and  real -coded  GA  up  to  sixth-order 
motion  (i.e.,  search  for  5  parameters)  and  plot 
the  results  in  Fig.  4(a).  It  is  observed  that  both 
binary-coded  GA  and  real-coded  GA-  have  much  . 
lower  complexity  than  exhaustive  search.  The 
difference  in  complexity  between  the  two  GAs  is 
only  minor.  The  projection  values  from  binary-coded 
GA  and  real-coded  GA  are  next.plotted  in  Fig.  4(b) 
as  circles  and  crosses,  respectively.  We.  see  that  the 
real-coded  GA  results  are  closer  to  die  truth  value 


of  1.0  than  the  binary-coded  GA  results,  especially 
for  higher  order  motions.  Since  binary-coded  GA 
searches  on  a  finite  grid  as  in  exhaustive  search,  it 
can  never  get  solutions  that  surpass  the  exhaustive 
search  result.  On  the  other  hand,  real-coded  GA  has 
the  ability  to  search  for  any  real  values  within  the 
search  range.  Consequently,  real-coded  GA  has  more 
chance  of  finding  a  better  solution.  The  same  trend  is 
also  observed  with  measurement  data  and  is  discussed 
further  in  the  next  section. 

V.  MEASUREMENT  DATA  PROCESSING 

We  next  apply  GA  on  some  measurement  data. 

The  data  were  collected  from  an  in-flight  aircraft 
using  ground  radar.  128  pulses  are  processed  to  form 
an  ISAR  image.  This  corresponds  to  an  imaging 
interval  of  about  2.5  seconds.  GA  is  evaluated  for  fine 
motion  compensation.  For  the  first  data  set,  the  image 
without  any  fine  motion  compensation  is  shown  in 
Fig.  5(a).  This  image  is  unfocused  due  to  the  residual 
motion.  We  select  a  range  cell  with  a  dominant 
scatterer  (range  cell  79)  and  apply  GA  to  estimate 
the  phase  based  on  a  third-order  translation  motion 
model.  The  resulting  images  from  binary-coded  and 
real-coded  GA  are  shown  in  Figs.  5(b)  and  5(c), 
respectively.  We  observe  that  the  two  GA  images 
are  as  focused  as  the  reference  image  in  Fig.  5(d) 
obtained  using  exhaustive  search.  The  corresponding 
projection  values  are  2401  and  2595,  as  compared 
to  the  exhaustive  search  result  of  2401.  We  continue 
this  comparison  using  19  other  imaging  intervals. 

Fig.  6(a)  shows  the  projection  values  (normalized 
with  respect  to  the  exhaustive  search  result)  for  the  20 
frames.  For  19  out  of  20  frames,  real-coded  GA  gives 
larger  projection  values  than  exhaustive  search.  The 
resulting  images  are  either  on  par  or  slightly  ,  better 
focused  than  those  obtained  from  exhaustive  search. 
For  binary-coded  GA,  7  frames  have  lower  projection 
values  and  are  of  inferior  image  quality  to  those 
from  exhaustive  search.  The  other  13  frames  have 
exactly  the  same  projection  values  as  the  exhaustive 
search  result  Similar  to  our  conclusion  earlier  based 
on  the  simulation  data,  our  experience  with  the 
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Fig.  5. 


(C)  (d)  •  .  r  - 

Tnmslatton  ‘  motion  compeasatkm  applied  to  measurement  data  using  a  third-order  motion  model,  (a)  Image  before  fee  motion 
compensation,  (b)  Binary-coded  GA  result,  (c)  Real-coded  GA  result.  (d)  Exhaustive  search  result 

•  J"”j 

.  .  •  *  .  /  A*  ' 


Fig.  6.  Comparison  of  performance  of  GA  and  exhaustive  search  using  measurement  data,  (a)  Projection  value,  (b)  Cdmputatico'  time: 


measurement  data  indicates  that  real-coded  GA 
consistently  outperforms  binary-coded  GA  in  terms 
of  accuracy.  The  computation  time  using  Matlab  codes 
on  a  Pentium  HI  750MHz  PC  is  summarized  in  Fig. 
6(b).  While  there  is  little  change  in  the  computational 
time  for  exhaustive  search  from  one  frame  to  another, 
the  times  for  binary-coded  and  real-coded  GA  exhibit 
large  Variations  in  these  single  run  results  due  to 
the  statistical  nature  of  GA.  The  average  times  for 
the  binary-coded  and  real-coded  GA  are  19.5  s  and 
11.5  s,  respectively.  This  is  compared  with  die  45.5 
s  from  exhaustive  search.  Finally,  we  note  that  for 
#  these  20  tones  the  third-order  model  is  adequate 


to  model  the  translation  motion.  Inclusion  of  higher 
order  translation  motion  or  rotational  motion  does  not 
improve  the  image  quality  for  this  data  set.  * 

For  a  second  data  set,  we  first  apply  third-order  1  * 
translation  motion  compensation  using  real-coded 
GA.  The  resulting  image  is  shown  in  Fig. '7(a).  It 
is  seen  that  the  selected  dominant  point  scatterer  at  ’ 
range  cell  64  is  not  well' focused.  This  means  that  • 
die  target  contains  higher  motion  that  cannot  be 
fully  compensated  by  the  third  order  motion  model. 
This  fact  is  also  revealed -by  the  spectrogram  of  the 
compensated  signal  in  Fig.  7(b),  as  .we  observe  a 
curved  joint  timeTrequency  (JTF)  trajectory  for 
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He.  7.  Motion  compensation  example  for  another  measurement  data  (a)  Image  after  third-order  translation  motion  compensation 
using  real-coded  GA.  (b)  Spectrogram  of  signal  in  range  cell  64. 


dwell  time 


.  .  •'»  >  ;  ..  « 

Pig.  8.  Higher-order  translation  motion  compensation,  (a)  Image  after  fourth-order  translation  motion  compensation  using  realTCOded 
GA.  (b)  Spectrogram  of  signal  in  range  cell  64.  (c)  Spectrogram  of  signal  in  range  cd!  71. 


the  scatterer.  Next,  a  fourth-order  motion  model  is  which  shows*  that  the  signal  trajectory  is  curyed  m  the 

tried  and  the  real-coded  GA  result  is  shown  in  Fig.  spectrogram.  Thus  rotational  motion  must  exist  in  this 

8(a).  From  this  figure,  the  reference  point  scatterer  data:  We  next  apply  fourth-order  rotational  motion  M 

is  better  focused  and  the  spectrogram  of  the  signal  compensation  using  real-coded  GA.  As  shown  in 

is  straightened  in  Fig.  8(b).  Fifth-order  translation  Fig.  9(a),  the  whole  target  is  much  better  focused  after 

motion  is*  also  attempted,  and  the  result  does  not  show  the  compensation.  The  spectrograms  of  the  signal  at 
much  improvement  However,  we  observe  that  other  both  range  cells  become  straightened  in  Figs.  9(b) 

point  scattered  in  Fig.  8(a)  are  still,  smeared.  This  and  9(c).  . While  real-coded  GA  takes  45  s  for  the 

is  confirmed  by  the  spectrogram  of  the  signal  at  a  phase  estimation  problem,  the  computation  time  for 

different  range  cell  (number.  71)  shown  in  Fig.  8(c),  fourth-order  exhaustive  search  is  estimated  to  be  over 
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fie.  9.  Higher-order  rotatkwul  motion  compcnsilkm.  (a)  Image  after  fbortb-anter  rotational  motion  compensation  using  real-coded  GA. 
(b)  Spectrogram  of  signal  in  range  cell  64.  (c)  Spectrogram  of  signal  in  range  cell  71. 


50  min  based  on  the  complexity  curve  in  Fig.  4(a). 
Therefore,  the  time  savings  of  GA  over  exhaustive 
search  is  quite  significant  in  this  real-world 
example. 

VI.  CONCLUSIONS 

In  this  paper,  GAs  have  been  applied  to  ISAR 
motion  compensation.  Based  on  the  AJTF  analysis, 

GA  is  used  in  the  phase  estimation  of  prominent  point 
scattered  on  the  target.  The  resulting  parameterized 
phases  are  then  used  for  translation  and  rotational 
motion  compensation.  Both  binary-coded  GA  and 
real-coded  GA  have  been  implemented  and  tested 
using  simulation  and  measurement  data.  It  is  found 
that  real-coded  GA  outperforms  binary-coded  GA  in 
terms  of  accuracy  in  the  phase  estimation  problem.  It 
is  also  shown  that  the  computational  complexity  of  the 
GA  search  is  much  less  than  that  of  exhaustive  search. 
The  time  savings  can  become  especially  significant 
when  the  target  exhibits  higher  order  motions. 
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A  new  technique  is  developed  for  range  alignment  In  inverse 
synthetic  aperture  radar  (ISAR)  Imaging.  The  shifts  made  to  the 
echoes  are  modeled  us  a  polynomial,  and  the  coefficients  of  this 
polynomial  are  chosen  to  optimize  a  global  quality  measure  of 
range  atignmenL  This  technkjoe  is  robust  agahkst  noise  and  target 
scintillation,  and  avoids  error  accumulation.  In  addition,  the  shift 
hi  the  time  domain  is  implemented  by  introducing  a  phase  ramp 
ia  the  frequency  domain,  which  removes  the  timitatico  of  Integer 
steps. 

L  INTRODUCTION 

Inverse  synthetic  aperture  radar  (ISAR)  utilizes  the 
Fourier  transform  to  resolve  the  scatterers  in  azimuth. 
Before  taking  the  Fourier  transform,  translation 
compensation  is  needed  to  remove  the  effects  of  the 
translation  between  the  radar  and  the  target  in  range. 
Translation  compensation  consists  of  range  alignment, 
which  shifts  the  echoes  such  that  the  signals  from 
the  same  scatterer  are  centered  at  the  same  range  ^ 
bin  in  different  echoes,  and  phase  adjustment,  which 
removes  die  Doppler  phase  caused  by  the  translation. 

If  no  prior  knowledge  is  available  about  the 
translation,  range  alignment  is  usually  based  on  the 
similarity  of  the  envelopes  of  the  echoes.  Typical 
methods  include  the  peak  method  [1],  the 
maximum-correlation  method  [1],  the  frequency- 
domain  method  [1],  die  Hough-transform  method  [2] 
and  the  minimum-entropy  method  [3], 

These  methods,  however,  have  a  variety  of 
disadvantages.  The  maximum-correlation  method,  for 
example,  aligns  each  echo  using  die  principle  that 
the  envelope  correlation  of  two  adjacent  echoes  is 
a  maximum  when  they  are  aligned.  Although  more 
robust  than  die  peak  method,  the  method  is  still 
somewhat  sensitive  to  noise  and  target  scintillation. 
Also,  it  has  the  defect  of  error  accumulation.  In 
addition,  the  shift  in  the  time  domain  has  the 
limitation  of  integer  steps,  which  means  that  even  if  . 
two  echoes  are  correctly  aligned,  there  may  still  be  an 
error  of  up  to  half  a  range  bin. 

We  develop  a  new  technique  for  range  alignment 
in  ISAR  imaging.  The  shifts  made  to  the  echoes  are 
modeled  as  a  polynomial,  and  the  coefficients  of  this 
polynomial  are  chosen  to  optimize  a  global  quality 
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the  proposed  codec  averagely  saves  8%  bit  rate  compared  to  JM3.9  with 
the  best  settings.  Complexity  comparison  between  JM3.9  with  five 
reference  frames  and  the  proposed  algorithm  is  also  analysed  The 
sprite  buffer  can  be  restricted  within  the  constant  size  (e.g.  2.25  times 
the  frame  size).  Therefore  fire  memory  cost  in  the  proposed  sprite  coding 
is  less  than  that  in  JM3.9.  For  computing  complexity  analysis,  we  only 
consider  the  computing  time  of  motion  estimation-  In  JM3.9,  the  local 
morion  estimation  (LME)  is  performed  five  times  (i.c.  once  for  each 
reference  frame).  In  the  proposed  codec,  LME  is  performed  only  once. 
The  total  rime  of  motion  estimation  in  the  proposed  codec  is  for  one 
LME  and  one  GME.  By  utilising  the  fast  GME  algorithm,  the  total  time 
for  LME  plus  GME  is  less  than  five  times  LME  in  JM3.9.  Similar  to  the 
traditional  dynamic  sprite  coding  techniques,  the  proposed  codec  has  the 
disadvantage  that  sprite  warping  has  to  been  performed  in  the  decoder. 
However,  considering  the  significant  improvement  of  coding  efficiency, 
the  extra  computing  complexity  for  sprite  warping  is  acceptable. 


Conclusions:  We  have  presented  a  highly  efficient  algorithm  for 
dynamic  sprite  coding.  The  high  coding  efficiency  is  achieved  due 
to  two  reasons.  First,  the  new  techniques  developed  in  JVT  codec  arc 
utilised;  secondly,  the  fractional  resolution  sprite  prediction  is  incor¬ 
porated  into  the  proposed  algorithm. 
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Shape  inversion  of  metallic  cavities  using 
hybrid  genetic  algorithm  combined 
with  tabu  list 

Yong  Zhou,  Junfei  Li  and  Hao  Ling 

An  Approach  combining  the  hybrid  genetic  algorithm  (GA)  with  the 
tabu  list  concept  is  proposed  to  iacrease  the  search  efficiency  of  the 
hybrid  GA  The  algorithm  is  applied  to  reconstruct  the  shape  of  a 
metallic  cavity  based  on  the  Ipswich  measurement  data.  Inversion 
results  show  good  agreement  with  the  actual  shape  and  significant 
improvement  m  convergence  rate  over  both  simple  GA  and  hybrid  GA. 


Introduction:  Electromagnetic  inverse  scattering  entails  five  recon¬ 
struction  of  the  shape  or  material  of  an  object  from  its  scattered  field 
data.  The  inverse  problem  can  be  cast  into  an  optimisation  problem 
whereby  the  difference  between  the  measured  fields  and  the  computed 
fields  from  a  forward  electromagnetic  solver  is  minimised.  Genetic 
algorithms  (GAs)  have  been  tried  as  the  global  optimisef  in  these 
problems  (1-3].  While  GA  is  well  suited  in  searching  for  the  global 
optimum,  it  suffers  from  slow  convergence.  Since  the  evolutionary 
process  for  the  standard  GA  to  reach  a  cost  minimum  is  in  general 
very  slow  in  comparison  to  a  local  search  algorithm,  a  natural 
improvement  to  speed  up  the  simple  GA  is  to  hybridise  the  simple 
GA  with  a  local  search.  This  type  of  algorithm  is  usually  called  the 
hybrid  GA  (HGA)  and  has  been  explored  by  researchers  i a  different 
disciplines  [4,  5].  While  showing  improvements  in  performance,  the 
hybridisation  of  the  two  approaches  also  causes  some  inefficiency.  As 
five  parent  selection  scheme  of  GA  gives  priority  to  the  best  members, 
it  usually  leads  to  a  population  that  is  highly  clustered  around  the 
local  minima.  This  clustering  is  necessary  for  tbo  simple  GA  to  evolve 
closer  to  the  exact  minimum.  For  HGA,  however,  since  the  local 
minima  have  been  completely  explored  by  the  local  search,  such 
clustering  will  lead  to  the  re-exploration  of  these  regions,  which  is 
quite  wasteful. 

Tabu  search  (TS)  is  another  global  search  strategy  that  has  been 
developed  for  combinatorial  problems  [6,  7].  It  is  a  local  search 
algorithm  with  memory.  The  most  important  feature  of  TS  is  that  it 
utilises  a  tabu  list  to  prevent  five  revisiting  of  local  minima.  In  this 
Letter,  we  propose  a  technique  combining  HGA  with  the  tabu  list 
concept  to  increase  the  efficiency  of  the  HGA.  The  tabu  list  is  adopted 
to  exclude  those  regions  in  the  parameter  space  that  have  already  been 
explored  by  the  local  search.  In  this  maimer,  there  will  be  no  revisiting 
of  the  explored  regions  and  the  GA  population  can  be  spread  out  to 
explore  new  regions,  finis  improving  the  search  efficiency.  We  apply 
this  algorithm  to  the  electromagnetic  inverse  problem  of  sliape 
reconstruction  of  metallic  cavity  structures  containing  strong  multiple 
scattering  effects.  Results  based  on  the  Ipswich  measurement  data  set 
are  presented. 


tabu  region 


parameter  space 


Fig.  1  Establishment  of  tabu  region  ’ 

HGA-tabu  approach:  In  our  HGA-tabu  implementation,  the  initial 
generation  is  produced  randomly.  The  new  population  is  then 
produced  through  the  selection,  crossover  and  mutation  operators. 
After  these  standard  GA  processes,  the  best  member  P  is  selected  as 
the  initial  guess  to  carry  out  a  local  search.  We  adopt  the  gradient 
search  reported  in  [8]  as  the  local  search  algorithm.  The  resulting 
local  minimum  in  the  parameter  space  is  denoted  as  Pmm  (sec  Fig.  1). 

is  then  placed  into  the  new  GA  population.  In  addition,  a 
gradient  search  is  also  carried  out  to  obtain  the  local  maximum 
from  the  same  initial  guess  m  order  to  estimate  the  extent  of 
the  local  minimum.  Once  both  local  searches  are  completed,  we  define 
the  region  that  is  centred  at  the  minimum  and  limited  by  the  radius 
as  the  tiabu  region’,  and  record  it  into  a  tabu  list. 
Symmetry  around  the  local  minimum  is  assumed  in  this  construct.  In 
subsequent  GA  reproductions,  all  of  the  new  members  are  checked 
against  this  tabu  list  to  ensure  that  none  is  in  the  tabu  regions  of  the 
sample  space.  Thus  file  population  is  forced  to  spread  out  to  the 
unexplored  regions,  resulting  in  higher  HGA  search  efficiency. 
Further,  a  new  tabu  region  is  appended  to  the  tabu  list  every  time  a 
new  local  minimum  is  explored  by  local  search. 

In  implementing  the  inverse  problem  to  reconstruct  the  shape  of  a 
cavity  from  its  scattered  field  data,  we  start  from  a  set  of  randomly 
created  shapes  that  are  described  by  N  ordered  points  in  a  two- 
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dimensional  apace.  The  profile  of  the  object  is  then  obtained  using 
spline  interpolation.  Next,  the  method  of  moments  (MoM)  solution  to 
the  electric  field  integral  equation  is  used  as  the  forward  electromag¬ 
netic  solver  to  generate  the  computed  scattered  field  Ecaf  from  each 
assumed  shape.  A  cost  function  is  defined  as  the  root-mean-squared 
(rms)  difference  between  and  the  measured  scattered-field  £*"“ 
The  HGA-tabu  algorithm  is  then  applied  as  the  optiraiser  to  minimise 
the  cost  function.  Binary-encoded  GA  is  used  in  our  implementation. 

Results'  We  have  applied  the  HGA-tabu  algorithm  to  reconstruct  the 
shape  of  a  metallic,  partially  open,  circular  cylindrical  cavity  with  a 
diameter  of  10.8  cm  (IpsOll  in  the  Ipswich  data  set)  [9j.  The 
measurement  was  taken  at  a  single  frequency  of  10  GHz  in  a  bistatic 
configuration.  There  were  a  total  of  36  transmitter  positions  around 
the  object  and  18  receiver  locations  for  each  transmitter  position.  The 
electric  field  was  parallel  to  the  axis  of  the  cylinder. 

The  number  of  the  population  for  GA  was  set  to  200,  the  geometry 
was  described  by  N—  5  points,  and  the  crossover  and  mutation  rates 
were  set  to  0.8  and  0.4,  respectively.  The  search  area  was  chosen  to  be 
16.2  x  162  cm.  We  first  tested  the  inversion  algorithms  using  MoM- 
simulated  field  data  as  the  input.  The  results  showed  that  the  HGA-tabu 
was  able  to  converge  to  the  correct  shape  after  an  average  of  75 
generations  and  the  final  shape  was  in  excellent  agreement  with  the 
actual  shape.  In  comparison  to  the  HGA,  the  HGA-tabu  also  showed  an 
improvement  of  about  100  generations  for  convergence. 


Fig.  2  Convergence  comparison  for  inversion  of  ipsOll  for  random 
search,  simple  GA.  HGA  and  HGA-tabu 


Fig.  3  IpsOll  inversion  results  from  measured  data 

a  Typical  taveiskni  result  by  smple  GA 
b  Typical  inversion  result  by  HGA-tabu 


Next,  we  applied  the  inversion  algorithms  to  the  actual  measured  data 
for  IpsOl  1.  Fig.  2  shows  the  convergence  comparison  between  random 
search,  ample  GA,  HGA  and  HGA-tabu.  All  the  results  were  averaged 
over  10  independent  runs  with  different  initial  populations.  As 
expected,  the  simple  GA  showed  improvement  over  the  random 
search.  The  HGA  further  improved  the  convergence  rate  of  the 
simple  GA-  The  best  results  were  consistently  obtained  by  the  HGA- 
tabu.  To  achieve  an  rms  of  0.55,  the  HGA  required  an  average  of  220 
generations  while  the  HGA-tabu  algorithm  required  only  an  average  of 
75  generations.  (We  note  here  that,  due  to  the  difference  between  the 
numerical  modelling  and  the  measurement,  the  rms  error  between  die 
MoM-computed  fields  from  the  exact  shape  and  the  measured  field  data 
is  0.73.)  Fig.  3a  shows  the  typical  shape  from  tbe  simple  GA  after  250 
generations  plotted  against  the  real  profile  of  the  cavity.  The  result 
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indicates  that  more  iterations  are  needed  for  convergence.  Fig.  3 b  shows 
tbe  typical  reconstructed  shape  from  the  HGA-tabu  after  75  genera¬ 
tions.  As  we  can  see,  the  inverted  shape  is  very  close  to  the  real  profile. 
The  overhead  of  implementing  the  gradient  search  in  each  generation  is 
about  10%  of  the  total  computation  cost.  The  time  for  the  tabu  list 
check  is  negligible,  as  there  is  no  cost  function  evaluation. 

Conclusion :  An  approach  combining  the  hybrid  genetic  algorithm 
with  the  tabu  list  concept  has  been  proposed  in  this  Letter.  The  tabu 
list  was  set  up  to  increase  the  search  efficiency  by  forbidding  revisits 
of  local  minima  already  explored  by  the  local  search.  The  algorithm 
has  been  applied  to  reconstruct  the  shape  of  a  metallic  cavity  based  on 
the  measured  Ipswich  data.  Invention  results  from  the  HGA-tabu 
showed  faster  convergence  and  higher  success  rate  than  those  of  the 
simple  GA  and  hybrid  GA.  The  computation  overhead  per  generation 
for  the  new  algorithm  was  small.  Tbe  algorithm  could  potentially  be 
useful  in  other  optimisation  problems. 
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Publicly  verifiable  authenticated  encryption 
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A  new  authenticated  encryption  scheme  with  public  verifiability  is 
presented.  The  new  scheme  requires  less  computational  costs  and 
communication  overhead  than  tbe  convention*!  signature- then- 
'  encryption  approaches.  Furthermore  the  message  is  not  divulged 
during  the  public  verification. 

Introduction:  Secure  and  authenticated  message  deliver/storage  is 
one  of  tbe  major  aims  of  computer  and  communication  security 
research.  Horster,  Michels  and  Petersen  (HMP  for  short)  [1]  proposed 
an  efficient  authenticated  encryption  scheme  with  lower  expansion 
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Sparse  parameterisation  of  electromagnetic 
scattering  data  using  genetic  algorithm 
with  adaptive  feeding 

J.  Li,  V.  Zhou  and  H.  Ling 

A  method  is  presented  to  parameterise  scattering  data  fora  complex 
targets.  Baaed  on  a  global  model  with  both  scattering  centres  and 
resonances,  a  genetic  algorithm  with  adaptive  feeding  is  proposed  fbr  a 
sparse  represenmbotL  The  algorithm  is  tested  with  measurcmail  data  and 
shows  better  performance  than  non -global  parameterisation  methods. 


process.  To  overcome  this,  a  parameterisation  based  on  a  community 
GA  with  adaptive  feeding  is  devised.  Fig.  1  illustrates  the  approach  in 
extracting  M  scattering  centres  and  N  resonances.  Each  solid  box 
represents  a  community  [10]  using  a  different  parameterisation  order 
number.  For  example,  the  highest  order  uses  M  scattering  centres  and  N 
resonances  while  the  lowest  order  uses  one  scattering  centre  and  no 
resonance.  The  parameterisation  consists  of  an  iteration  process  as 
■follows.  First,  for  each  community,  standard  GA  operations  including 
selection,  crossover  and  mutation  are  used  to  reproduce  members  in  the 
next  generation  for  better  solutions.  Secondly,  at  the  end  of  each 
generation,  the  residual  signal  of  each  community  is  calculated  as  die 
error  between  the  best  solution  in  die  community  and  the  original  data 
£*"(/)  and  is  parameterised  with  GA.  The  order  number  for  residual 
parameterisation  is  specified  in  the  dashed  box.  It  is  die  difference 
between  the  order  number  of  die  current  community  and  the  next  higher 
community.  Thirdly,  the  parameters  from  the  best  solution  of  a  lower 
order  community  and  its  residue  arc  combined  to  form  a  candidate 
solution  in  the  next  higher  order  community.  A  zero-mean  Gaussian 
perturbation  is  added  during  this  step  to  create  a  community-level 
mutation.  By  adaptively  feeding  the  solutions  from  the  lower  order 
communities  forward  to  the  higher  order  communities,  the  convergence 
of  the  highest  order  community  is  significantly  accelerated  without 
sacrificing  the  optimality  of  the  final  solution. 


Introduction:  Obtaining  a  sparse,  physical  representation  of  electro¬ 
magnetic  scattering  data  from  a  complex  target  is  a  problem  of 
fundamental  importance  in  radar  signature  analysis  [1-9].  The  scat¬ 
tering  centre  model  is  the  standard  way  to  represent  scattering  from 
large  targets  and  has  been  used  with  success  by  the  radar  signature 
community  for,  over  two  decades.  Many  techniques  including  super- 
resolution  [1,  2],  CLEAN  [3],  genetic  algorithms  (GA)  [4-b]  and 
evolutionary  programming-based  CLEAN  [7]  have  been  reported  for 
determining  model  parameters  based  on  the  scattering  centre  model. 

For  targets  containing  convex,  interior  structures  such  as  cavities,  a 
model  combining  scattering  centres  and  resonances  has  been  shown  to 
be  a  more  efficient  and  physically  meaningful  representation  of  both 
exterior  and  interior  scattering  features  [8,  9].  However,  finding  the 
model  parameters  in  such  cases  is  more  challenging,  since  the  scattering 
centre  and  resonance  bases  have  complementary  behaviours  in  time  and 
frequency.  In  [8],  a  CLEAN-bascd  algorithm  was  used  to  extract  one 
scattering  centre  and/or  resonance  at  a  time  iteratively.  In  [9],  Prony’s 
method  was  first  used  to  extract  all  the  scattering  centres  and  then  all 
the  resonances.  One  drawback  of  these  methods  is  that  the  parameter! - 
satioc  results  are  not  very  sparse  since  die  scattering  centres  and 
resonances  are  not  extracted  simultaneously. 

To  improve  the  sparsity,  wc  present  in  this  Letter  a  global  algorithm 
to  parameterise  complex  scattering  data  using  the  combined  scattering 
centre  and  resonance  model.  The  method  is  based  on  a  GA  with 
adaptive  feeding.  The  latter  is  devised  to  compensate  for  the  disparity 
in  strength  between  scattering  centres  and  resonances  and  improve  the 
performance  of  the  GA. 


GA  with  adaptive  jeeding:  The  scattering  model  is  assumed  to 
comprise  responses  from  both  scattering  centres  and  resonances  as 

[9]: 


E{f)  =  L  a~r'K'PllJ + 5  -/,) + a 


-ftrfx. 


0) 


where  M  and  N  arc  the  number  of  scattering  centres  and  resonances, 
respectively,  and  f  the  frequency.  For  each  scattering  centre,  im  is  the 
time  delay,  a*  the  frequency  dependency  coefficient  and  am  the 
complex  amplitude.  For  each  resonance  of  complex  amplitude  bm  fn 
is  the  resonant  frequency,  tm  the  tum-on  time  and  the  {2-factor.  The 
parameterisation  process,  can  be  formulated  as  a  minimisation  problem: 

/*,/«■  ft,,  tw}  =  argmin  UWO  -  E*if)h  (2) 

where  £*"  denotes  the  measurement  data  to  be  parameterised.  The 
amplitudes  am  and  bH  are  not  included  in  the  bracket  as  they  can  be 
derived  from  other  unknowns  using  minimum  least  squares  fitting. 

GA  has  been  used  in  many  engineering  applications  as  a  global 
optimisation  scheme.  However,  here  we  find  that  the  standard  simple 
GA  (SGA)  has  difficulty  m  converging  to  the  desired  global  optimum. 
Since  the  energy  in  a  resonant  term  is  typically  much  lower  than  that  in 
a  scatterer  centre,  the  resonant  terms  are  easily  missed  in  the  SGA 
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Fig.'  1  GA  with  adaptive  feeding 

Best  solution  from  lower  order  community  (solid  box)  and  residue  (dashed  box) 
combined  and  fed  into  next  higher-order  community.  Convergence  of  highest 
order  community  with  M  scattering  centres  and  N  resonances  accelerated 


Fig.  2  Comparison  of  three  pammeterisaiion  results  of  VFY  2J8  measure¬ 
ment  data 

- CLEAN  with  scattering  centres  only 

- CLEAN  with  scattering  centres  and  resonances 

- -  GA  with  adajrtive  feeding 

VYF  2J8  measurement  data  results:  The  algorithm  was  first  tested 
using  numerical  simulation  data  from  a  well-understood  target,  a  plate 
with  a  partially  open  cavity  [9].  The  proposed  method  successfully 
extracted  the  four  dominant  scattering  centres  and  three  resonances 
with  a  5%  RMS  error.  By  comparison,  the  RMS  error  from  the 
CLEAN  method  is  10%  after  20  terms,  while  the  SGA  always  missed 
the  weakest  resonance.  We  next  applied  it  to  the  VFY  218  measure¬ 
ment  data  [11].  The  scattering  data  came  from  a  1:30  scale  model 
aircraft  using  horizontal  polarisation  in  the  8  to  16  GHz  frequency 
band.  The  look  angle  was  at  19.6°  from  nose-on  so  that  the  inlet 
contribution  was  prominent  in  the  return.  Fig.  2  shows  RMS  error 
against  total  model  order  number  (M+N)  fbr  three  different  meth¬ 
ods.  The  GA  results  were  averaged  over  six  runs.  The  CLEAN  curve 
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with  scattering- centre- only  model  decreases  very  little  after  30  terms, 
indicating  this  model  is  very  inefficient  in  modelling  the  resonance 
part  of  the  data  after  the  scattering  centres  are  extracted,  The  CLEAN 
curve  with  both  scattering  centre  and  resonance  model  is  better. 
However,  the  rate  of  convergence  still  slows  down  considerably 
after  the’  first  18  terms.  The  GA  curve  shows  the  best  sparsity.  It 
requires  only  20  terms  to  achieve  the  same  accuracy  as  the  CLEAN 
approach  with  40  terms.  The  model  orders  used  in  the  GA  are  M  *  14 
and  V=6.  However,  we  found  that  the  results  were  not  very  sensitive 
to  the  model  order  selection. 


frequency,  GHz 
b 

Fig.  3  Accuracy  of  parameterization  in  time  and  frequency  domains 

a  Tune  domain 
b  Frequency  domain 


To  further  interpret  the  physical  significance  of  the  GA-pani- 
metcrised  data,  we  correlated  die  extracted  scattering  centre  positions 
with  the  peaks  in  the  target  range  profile  and  found  that  they  lined  up 
well.  Furthermore,  the  two  strongest  resonances  extracted  are  at 
frequencies  of  9.8  and  1 1.3  GHz.  This  is  consistent  with  die  size  of 
die  rectangular  engine  inlet  openings,  which  have  dimensions  of 

2.5x15  cm.  (The  cutoff  frequencies  of  the  TEoi  and  TEu  modes 
are  estimated  at  10  and  11.7  GHz,  respectively.)  The  other  four 
resonances  at  8.6,  9.1,  9.4  and  13.3  GHz  are  harder  to  interpret  given 

die  complex  shape  of  the  actual  inlet  structure. 

Comparisons  of  the  parameterised  result  with  the  original  measure¬ 
ment  data  in  the  time  and  frequency  domains  are  shown  in  Fig.  3.  We 
see  fairly  good  agreements  between  the  two.  We  suspect  the  small 
parameterisation  error  to  be  due  to  die  model  mismatch  of  (1)  to  the 
complex  measurement  data.  Thus,  increasing  the  model  order  for  this 
data  does  not  reduce  the  error  significantly.  We  also  processed  data  from 
0°  to  180°  from  nose-on  in  5°  increments,  and  found  foe  GA  with 
adaptive  feeding  to  consistently  outperform  CLEAN  at  all  angles. 


Conclusions:  We  have  proposed  a  GA-based  method  for  parameters 
sing  scattering  data  from  complex  targets.  Based  on  a  global  model 
with  both  scattering  centres  and  resonances,  our  method  uses  GA  with 
the  adaptive  feeding  idea  to  simultaneously  extract  all  the  model 
parameters.  The  proposed  method  can  achieve  sparser  results  than 
other  non-global  based  methods.  The  effectiveness  of  the  proposed 
method  is  demonstrated  using  the  VFY  218  measurement  data.  The 
resulting  sparse  model  facilitates  target  feature  interpretation  and  can 
be  used  for  signature  reconstruction  in  modelling  and  simulation 
applications. 
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400  mW  uncooled  MiniDIL  pump  modules 

S.  Mohrdiek,  T.  Pliska,  R.  Battig,  N.  Matuschek, 

B.  Valk,  J.  Tioger,  P.  Mauron,  B.E.  Schmidt, 

LD.  Jung,  C.S.  Harder  and  S.  Enochs 

A  new  generation  of  wavelength  stabilised,  uncookd  980  run  pump 
modules  in  MiniDIL  housings  is  presented,  enabling  400  mW  ex-fibre 
power  over  a  temperature  range  of  10°C  to  70°C.  At  100°C  200  mW 
power  is  still  obtained  with  a  robust  fibre  coupling  scheme. 

Introduction:  As  the  focus  in  optical  telecommunications  systems 
turns  more  towards  affordability,  there  is  a  push  to  produce  EDFAs  of 
lower  cost,  smaller  size  and  less  power  consumption.  Operation  of 
980  tun  pump  modules  without  a  thermo-electric  Cooler  (TEC)  has 
been  presented  in  (1].  The  removal  of  the  bulky  and  power-consuming 
TEC  allowed  us  to  develop  pump  modules  in  a  smaller,  less  expensive 
MiniDIL  housing.  Though  low  cost  is  crucial,  performance  and 
reliability  comparable  to  conventional  Butterfly-type  modules  has  to 
be  demonstrated,  in  order  to  satisfy  foe  yet  stringent  requirements  for 
metro  systems. 

In  this  Letter  we  present  results  of  550  mW  fibre  coupled  power  at 
25°C,  400  mW  at  70°C  and  200  mW  at  the  extreme  temperature  of 
100°C,  with  Minimi,  modules  incorporating  the  latest  developments  in 
pump  laser  devices  [2,  3]  and  wavelength  stabilisation  by  fibre  Bragg 
gratings  (FBGs)  [4].  Little  change  in  fibre  coupling  efficiency  with 
temperature  demonstrates  the  robustness  of  the  fibre  alignment  scheme. 
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Application  of  adaptive  chirplet  representation  for 
ISAR  feature  extraction  from  targets  with  rotating 

parts 


J.  Li  and  H.  Ling 


Abstract:  The  problem  of  feature  extraction  from  inverse  synthetic  aperture  radar  (ISAR)  data 
collected  from  targets  with  rotating  parts  is  addressed.  In  traditional  ISAR  imaging,  rigid-  y 
motion  is  usually  assumed.  When  non-rigid-body  motions  are  present,  it  is  not  possible  to  obtain  a 
focused  image  of  both  the  target  and  the  rotating  pan.  To  solve  this  problem  the  radar  signal  is  hrs 
parameterised  using  the  adaptive  chirplet  signal  representation.  The  signal  from  the  body  and  that 
from  the  rotating  pan  are  then  separated  in  the  parameter  space.  Point-scatterer  simulation  results 
show  that  belter  geometrical  features  of  the  body  and  better  nucro-Doppler  features  ol  the  rotating 
pan  can  be  extracted  after  the  separation.  The  algorithm  is  also  demonstrated  using  the 
measurement  data  from  an  in-flight  aircraft  and  a  walking  person. 


1  Introduction 

Recently,  there  has  been  increasing  interest  in  studying  the 
so-called  micro-Doppler  phenomenon  [1,2]  for  radar  target 
identification  applications.  Micro-Doppler  is  used  to 
describe  the  fine  Doppler  feature  from  some  moving  part 
on  the  target  that  is  different  from  the  main  body  Doppler 
feature.  In  most  of  the  conventional  work  on  inverse 
synthetic  aperture  radar  (ISAR)  imaging,  the  target  is 
assumed  to  have  rigid-body  motion  [3,  4].  However,  non- 
ri°id-body  targets  can  often  be  found  in  real-world 
stations.  As  a  simple  case,  a  target  may  consist  of  a 
main  body  and  a  rotating  part.  For  example,  an  in-flight 
aircraft  with  jet  engine  rotation,  a  ship  with  scanning 
antenna  motion  and  a  ground  vehicle  with  spinning  tyre 
motion  all  involve  this  type  of  configuration.  Under  these 
conditions,  difficulties  in  understanding  the  resulting  ISAR 
image  arise  due  to  the  violation  of  the  rigid-body 
assumption. 

In  this  paper,  we  set  out  to  extract  better  target  features 
from  ISAR  data  when  a  target  has  a  rotating  part  beside  the 
main  body.  The  challenge  is  that  the  body  image  is 
contaminated  due  to  the  interference  irom  the  rotating 
part.  It  is  also  more  difficult  to  extract  the  motion 
information  from  the  rotating  part  as  it  is  overshadowed 
by  the  body  returns.  Our  approach  is  to  first  parameterise  the 
radar  signal  using  the  adaptive  chirplet  representation  [5, 6], 
The  chirplet  basis  is  a  four-parameter  function  localised  in 
the  joint  time- frequency  plane.  While  both  Gaussian  [7] 
and  chirp-type  [8]  bases  have  been  reported  for  joint 
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time -frequency  processing  of  ISAR  data,  the  chirplet  basis 
is  selected  to  represent  the  radar  signal  in  this  paper.  Since 
both  amplitude  modulation  (AM)  and  frequency  modulation 
(FM)  are  part  of  the  basis,  the  chirplet  can  more  efficiently 
represent  the  radar  signal  from  a  target  with  a  rotating  part. 
With  the  adaptive  chirplet  representation,  different  motion 
behaviours  of  the  target  components  are  mapped  into 
different  parameters  of  the  corresponding  bases.  Conse¬ 
quently,  the  returns  from  the  body  and  the  rotating  part  can 
be  more  easily  separated.  After  the  separation,  better  target 
feature  extraction  can  be  realised  by  processing  the  two 
parts  individually.  This  includes  both  the  extraction  of  the 
geometrical  features  from  the  main  body  and  the  micro- 
Doppler  features  from  the  moving  part. 

In  the  following  section,  we  present  the  model  and 
formulation  of  the  problem.  After  a  close  examination  of  the 
point-scatterer  signal  model  of  a  target  with  individual 
motions,  we  show  that  the  chirplet  basis  is  well  suited  for 
parameterising  and  separating  the  rotating  part  signal  from 
the  main  body  signal.  The  chirplet-based  adaptive  signal 
representation  algorithm  is  tested  with  point-scatterer 
simulation  data  and  results  are  shown  from  two  sets  ol 
measurement  data.  The  first  data  set  is  from  an  in-flight 
aircraft  with  jet  engine  rotation  motion.  The  second  data  set 
is  from  a  walking  person  with  arm  swinging  motion. 

2  Signal  model  and  formulation 

2. 1  Point-scatterer  model  of  radar  signal  from 
target  with  rotating  part 

The  point-scatterer  model  is  usually  used  in  radar  imaging 
to  model  the  radar  signal  scattered  by  an  unknown  target. 
In  this  model,  the  radar  return  signal  is  expressed  as  a  sum 
of  point-scatterer  responses 

£(/,  /)  =  om  exp  |  -j~r  [Ain (0  +  x»>  cos  M 

+y,,,sin0m(/)]j  (1) 

/££  Proc.-Ritdctr  Sonar  Navig..  Vol.  150 ,  No.  4.  August  2003 


where  the  radar  signal  £  is  a  two-dimensional  function  of 
transmitting  radar  frequency  /  and  pulse  dwell  time  /. 
The  target  consists  of  M  poinl-scatterers,  each  with  position 
(.v„„  y,„)  and  complex  scattering  coefficient  anr  Suppose  the 
radar  is  stationary,  the  target  motion  is  described  by  the 
translation  motion  Rm(f)  and  the  angular  motion  for 
each  scatterer. 

A  rigid-body  target  is  usually  assumed  in  traditional 
1SAR  imaging,  i.e.  all  the  point -sc atterers  in  (1)  share  the 
same  translation  motion  /?,„(/)  and  rotational  motion  9m(t). 
Here,  we  shall  consider  a  non-rigid-body  target  consisting  of 
two  parts,  a  main  body  and  a  rotating  part.  In  this  case,  we 
can  simplify  the  model  in  (1)  by  using  different  motions  for 
the  two  parts  while  still  applying  the  rigid-body  assumption 
lor  each  part.  This  leads  to 

£(/,  /)  =  £/?(/,  /)  +  £/?(/ 1 0 

=  ^2  a m  exp  |-7— [£/i(0  +  xm  cos  0B(t) 

m=l  ^  C 

}n  f  A'jtf 

H=\  ' 

+.V„  cos  4(?)  +  y„  sin  4(f)]  |  (2) 

with  subscripts  B  and  R  denoting  the  body  and  the  rotating 
part,  respectively. 

Both  the  main  body  and  the  rotating  part  move  with 
respect  to  the  radar.  The  difference  is  that  the  rotating  part 
has  an  additional  rotation  motion  beside  all  the  motions  of 
the  main  body.  For  the  main  body,  during  the  imaging 
interval  we  can  apply  the  small-angle  approximation  usually 
used  in  ISAR  imaging.  That  is 

f  cos  0B(i)  «  1 

\  sin  4(0  M  4(0  (3) 

We  also  assume  that  a  standard  motion  compensation 
algorithm  [4,  9,  10]  has  been  utilised  to  remove  both  the 
translation  motion  and  the  nonuniform  rotational  motion 
from  the  body,  after  which  we  can  write 

\  0B{t)  <-  toBi  (4) 

where  (oB  is  the  effective  body  rotation  rate  after  the  motion 
compensation.  The  arrow  symbol  above  is  used  to  indicate  a 
new  assignment  of  the  variable  on  the  left  after  the  motion 
compensation  operation. 

For  the  rotating  part,  the  motion  relative  to  the  main  body 
is  rotation  only.  This  implies  that  the  rotating  part  has  the 
same  translation  motion  as  that  of  the  body  while  the 
rotation  motion  will  change  accordingly,  i.e. 


However,  the  small-angle  approximation  does  not  hold  for 
the  rotating  part.  Since  its  rotation  rate  is  usually  much 
larger  than  that  of  the  main  body,  a  rotating  scatterer  might 
undergo  many  cycles  while  the  main  body  rotates  only  a  few 
degrees  during  the  imaging  interval.  Substituting  (3)-(5) 
into  (2),  we  have 

m  r  4  nf  } 

£X/v)=X/T"'expi  ~j~r  tv"'+ymC0B,i  r 

m— i  ^  * 

+y^(j„exp  | -i~~r [*v» cos^(0 +v ,, sin  0*  ( / )]  j  (6) 


which  is  the  radar  signal  from  a  target  with  a  rotating  part  in 
the  two-dimensional  (frequency,  dwell  time)  domain  after 
motion  compensation  of  the  main  body.  Since  it  is  more 
efficient  to  process  range  compressed  data,  we  Fourier 
transform  (6)  with  respect  to  / and  bring  the  radar  data  into 
the  (range,  dwell  time)  domain.  The  radar  signal  through  a 
fixed  range  cell  r  is  given  by 


Er{t)  =  y  a«>  4»sinc  I  (r  -  -v»>)  I 

X  exp  |  +j^r  ('■  -  xm  -  .W)} 

+  Y  onfb* sincj^— ^  [/•  -  .Y„  cos  4(0 
/;=1  ^ 

-  v„  sin  4(f)] |  exP  |  +j~~ lr  -  x»  cos  ^*(0 
-y„  sin  4(01 1 

where/;.  and/},,,.  are  the  carrier  frequency  and  the  bandwidth 
of  the  radar,  respectively. 

Some  observations  can  be  made  here  about  (7).  There 
exist  substantial  differences  between  the  main  body  signal 
and  the  rotating  part  signal.  Each  body -scatterer  in  the  first 
term  has  constant  amplitude  am  and  constant  Doppler 
frequency  -(2 /i/c)(%yw  with  respect  to  /.  However,  the 
signal  of  each  rotating  scatterer  in  the  second  term  contains 
both  AM  and  FM  components.  This  can  be  seen  by  the 
presence  of  the  time- varying  function  0*(/)  in  both  the  sine 
and  the  exponential  terms.  Consequently,  a  second  Fourier 
transform  of  (7)  with  respect  to  t  will  focus  the  target  body 
in  cross-range,  but  not  the  rotating  part  return.  This  results  in 
the  observed  interference  from  the  rotating  part  in  the  ISAR 
image. 

2.2  Chirplet  basis 

The  chirplet  basis  function  [5,  6]  is  well  suited  for 
parameterising  the  AM-FM  radar  signal  in  (7).  A  chirplet 
is  a  four-parameter  basis  of  the  form 

\4 

hk  (?)  =  Q0  4  exp  {-aA  (?  -  ?*)2}exp  {-j2irfk(i  -  ?*) 

where  tk  is  the  time  centre  of  the  signal ,/*  is  the  centre 
frequency,  0k  is  the  frequency  modulation  rate  and  ak 
defines  the  time  extent  of  the  signal.  The  joint  time- 
frequency  plot  of  a  chirplet  function  is  illustrated  in  Fig.  la. 

Actually,  the  chirplet  basis  is  one  of  the  many  options  that 
can  be  used  to  model  the  radar  signal  accurately.  However, 
there  are  some  attractive  attributes  of  this  basis.  First,  the 
basis  function  is  an  AM-FM  signal  and  only  a  sparse  set  of 
these  bases  is  needed  to  approximate  the  time-frequency 
structure  of  the  radar  signal  in  (7).  Secondly,  the  chirplet 
basis  is  a  well  understood  basis  with  only  four  parameters. 
Only  moderate  computation  time  is  needed  to  search  for 
the  basis  parameters.  Thirdly  and  most  importantly,  the 
parameters  of  the  chirplet  can  be  used  to  separate  the  two 
components  of  the  signal.  This  is  because  signals  from  the 
main  body  and  the  rotating  part  are  captured  by  chirplet 
bases  with  different  parameters. 

To  see  this  more  explicitly,  let  us  assume  a  first-order 
rotational  motion  m  the  time  neighbourhood  of  each  chirplet 

0B(i)  =  Bc  4-  toR(t  -  tk)  (9) 


.w)} 


‘/i  cos  4(0 


-  V„  Sin  4 


,v„  cos  4(0 


l EE  Proc.-Rudar  Sonar  Navig..  Voi  150.  No.  4.  August  2003 


285 


Fig.  1  Cliirplel  basis  and  chirplet  parameters 

a  Joint  lime-frequency  representation  of  a  chirplet  basis 
/.  Distribution  of  chirplet  parameters  for  the  mam  body  (solid) 
parts  (doited)  and  separation  thresholds  (dashed) 


and  rota  lion 


where  0,  is  the  angle  at  the  time  centre.  The  rotating  part  is 
assumed  to  have  a  constant  rotation  rate  o>K  during  the  time 
interval  near  it  although  it  could  have  more  complex 
motions  during  the  whole  imaging  interval. 

After  substituting  (9)  into  (7).  we  take  the  first  and  second 
derivatives  of  the  phase  term  with  respect  to  /  and  compare 
them  to  those  from  (8)  to  arrive  at  expressions  for./*  and  ft. 
The  results  can  be  written  as 

f.  =  -  —  <»n[.v„  sin  <>jk(i  -  >k)  -  cos  wtti‘  ~  '*)! 
c 

=  -  —  ln<t)R  sin  [<»K(t  -  tk)  -  c»]  ( 

c 


j3k  =  -“7  ln(oRcos  K(/  -  tk)  -  C„]  (1W?) 

where  (/„,  l„)  are  the  polar  representations  of  (.v/n  v  ).  It  can 
be  seen  from  (10a)  and  (10/?)  that  the  parameters  /*  and  pk 
are  distributed  along  an  ellipse  as  follows: 


Jt  ,  & 

(Ih(»rY  \C  y 


00 


where  the  size  and  the  axial  ratio  of  the  ellipse  are  controlled 
by  i»R  and  /„.  Similarly,  the  equation  to  associate  the 
chirplet  parameters  with  the  main  body  signal  is  given  by 


(12) 


where  the  scatterer  radial  length  /,„  and  the  rotation  rate  coD 
are  used  for  the  main  body. 

Even  though  (11)  and  (12)  have  exactly  the  same  form, 
the  main  body  and  the  rotating  pan  are  separable  in  the 
parameter  space  because  of  their  diflerenl  motions. 
Essentially,  while  the  sizes  of  the  two  pans  are  comparable, 
the  rotating  part  rotates  much  faster  than  the  mam  body 
during  the  imaging  interv  al,  i.e. 

(0R  >  ('?//  O3) 


Consequently,  the  chirplet  parameters  /*  and/?*  for  the  main 
body  and  the  rotating  part  are  distributed  very  differently  in 
the  parameter  space.  A  rotating  part  scatterer  is  represented 
as  a  larger  and  rounder  ellipse  while  a  body- scatterer  is 
represented  as  a  smaller  and  flatter  ellipse.  Actually,  the 
ellipse  for  the  main  body  is  nearly  a  line  segment  on  the 
fk  axis  since  the  first  term  of  (7)  is  assumed  to  have  zero 
Doppler  rate.  The  different  distributions  of  the  chirplet 
parameters  are  illustrated  in  Fig.  1/?,  where  the  outer  ellipse 
represents  the  rotating  part  signal  while  the  inner  one 
represents  the  main  body  signal.  A  simple  criterion  to 
separate  the  two  parts  can  thus  be  defined;  the  body  signal 
has  small./*  and  /?*  while  the  rotating  part  signal  has  either 
large./*  or  /?*. 

Another  interesting  observation  we  can  make  from  the 
above  discussion  is  that  the  main  body  and  the  rotating  part 
signals  have  large  overlaps  in  jk.  while  they  have  little 
overlap  in  ft.  Therefore,  the  Doppler  rate  is  more  important 
than  the  Doppler  frequency  in  separating  the  two  signal 
components.  This  point  will  be  further  illustrated  by 
examples  later. 

To  summarise,  if  we  parameterise  the  radar  signal  in 
question  into  a  set  of  chirplet  bases,  it  is  possible  to  separate 
the  contributions  from  the  target  body  and  the  rotating  part 
based  on  the  parameters  of  the  chirplet  bases,  as  we  have 
discussed  above. 


2.3  Signal  separation  based  on  adaptive 
chirplet  signal  representation 
To  decompose  the  radar  signal  into  a  set  of  chirplet  bases, 
we  apply  the  adaptive  signal  parameterisation  algorithm 
f  1 1  12].  We  start  with  the  radar  signal  in  a  fixed  range  cell 
with  returns  from  both  the  body  and  the  rotating  part,  which 
is  labelled  as  Er{t)  in  (7).  Next,  we  parameterise  Er(t)  by 
projecting  the  signal  onto  chirplet  bases  ol  diflerenl 
parameters  and  find  the  one  with  the  maximum  projection 
value.  Next,  a  residual  signal  is  generated  by  subtracting  the 
contribution  of  the  just-found  basis  from  the  signal.  This 
process  is  then  iterated  to  generate  a  series  of  chirplet  basis 
functions  that,  when  summed,  can  approximate  the  original 
signal.  The  steps  are  summarised  below: 

Step  1.  Set  iteration  index  number  k  to  1  and  the  residual 

signal  Rk(t)  to  E,-(t)  ... 

Step  2.  Find  the  Mh  chirplet  ltk(t)  by  maximising  the 
projection  from  the  residual  signal  Rk(t)  onto  the  basis,  i.e. 

{;*,/*,  a*, ft}  =  arg  max  \{Rk(t),hk(t))\  (14) 

where  the  inner  product  is  defined  as 

( Rk(i)M>))=  r  *»««(»)*  05) 

jin 

The  radar  data  are  assumed  to  exist  over  the  time  interval  t{) 
to  /,.  The  coefficient  of  the  chirplet  is  the  corresponding 
projection 

c*  =  {Rk(i)M>))  O6) 

Step  3.  Subtract  the  extracted  signal  from  the  residual 

Rk+l(t)*-Rk(<)-cM<)  (,7) 

Step  4.  Increment  k  by  one  and  repeat  steps  2  and  3  until  k 
reaches  a  preset  number  or  until  the  energy  of  the  residual 
signal  is  below  some  threshold  set  based  on  the  signal- 
to-noise  ratio.  Suppose  N  chirplets  are  found  from  this 
procedure,  the  radar  signal  is  parameterised  as 
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(18) 

After  the  parameterisation  of  the  radar  signal,  the  body 
signal  can  be  separated  from  the  rotating  part  signal  using 
the  criteria  discussed  previously.  We  classify  those  chirplets 
with  small/*  and  pk  as  the  main  body  components  and  the 
chirplets  with  either  large/-  or  large  //  as  the  rotating  part 
components.  The  final  body-only  signal  and  the  rotating  part 
signal  are  assembled  from  the  corresponding  chirplet  bases 
according  to  (18). 

Following  the  separation,  we  can  process  the  main  body 
signal  and  the  rotating  part  signal  individually  for  better 
information  extraction.  Based  on  our  discussion  about  (7), 
for  the  target  body  the  feature  of  interest  is  the  geometrical 
information  in  the  1SAR  image.  A  better  body  image  can  be 
reconstructed  after  removing  the  rotating  part  components. 
For  the  rotating  part  signal,  it  may  be  impossible  to  also 
construct  a  focused  image  of  the  rotating  part  if  the  PRF  of 
the  radar  is  too  low.  However,  it  is  possible  to  extract  useful 
information  about  the  motion  ol  the  rotating  part  front  the 
separated  data. 

3  Point-scatterer  simulation  results 

We  first  test  our  algorithm  with  point-scatterer  simulation 
data.  Six  point-scatterers  are  used  in  the  simulation  with  five 
points  representing  the  rigid  body  and  one  representing  the 
rotating  part.  The  positions  and  the  strengths  of  the  six 
scatterers  are  shown  in  Fig.  2.  Scatterer  6  rotates  around 
scatterer  2  at  a  rate  of  6.67  Hz  and  a  rotation  radius  of  20  cm. 
We  assume  the  radar  has  a  10  GHz  centre  frequency, 
800MHz  bandwidth  and  1400Hz  PRF.  The  target  body 
rotates  about  4°  over  384  pulses  during  the  data  collection 
time. 

Simulated  radar  data  are  generated  using  the  point- 
scatterer  model  in  (6).  The  resulting  radar  image  is  shown 
in  Fig.  3.  The  three  point-scatterers  in  the  centre  range 
cell  are  shadowed  by  a  noisy  vertical  micro-Doppler  band 
due  to  the  motion  of  the  rotating  point-scatterer.  Our 
objective  is  to  reconstruct  the  five  body-sealterers  and  to 
estimate  the  rotation  rate  of  the  rotating  scatterer  from  the 
radar  signal. 

Different  behaviours  of  the  body  and  the  rotating 
part  are  better  identified  in  the  joint  time- frequency 


X 


Fig.  2  Point-scatterer  representation  of  the  original  target 
consisting  of  five  rigid  points  (1—5)  and  one  rotating  point  (6) 
w  ith  strengths  2.  5.  2.  1,  L  and  3.33,  respectively 
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Fig.  3  Simulated  ISAR  image  of  the  target  with  non-rigid-body 
motion 


Fig.  4  Spectrogram  of  the  radar  signal  through  range  cell  65 


domain.  The  spectrogram  obtained  irom  the  short-time 
Fourier  transform  is  shown  in  Fig.  4  using  the  data  in  range 
cell  65,  which  contains  responses  from  scatterers  I,  2,  3 
and  6.  In  this  Figure,  we  see  interesting  features  about  the 
target.  First,  there  are  three  horizontal  Doppler  lines.  The 
one  at  zero  Doppler  is  due  to  scatterer  2.  The  two  at 
±100  Hz  are  due  to  scatterers  1  and  3.  Secondly,  there  is  a 
sinusoidal-like  micro-Doppler  curve  due  to  the  rotating 
scatterer  6.  Amplitude  modulation  of  this  signal  is  also 
observed. 

Following  the  steps  in  Section  2.3,  we  first  parameterise 
the  signal  using  N  -  100  chirplets.  The  spectrogram  of 
the  resulting  parameterised  signal  is  shown  in  Fig.  5.  We 
see  fairly  good  agreement  between  the  original  signal  and 
the  parameterised  signal.  Next,  we  separate  the  contri¬ 
butions  from  the  static  and  dynamic  parts  of  the  target 
based  on  the  Doppler  frequency  fk  and  the  Doppler  rate  pk 
of  the  chirplet  bases.  A  simple  threshold  of  3200  Hz/s  on 
the  Doppler  rate  and  300  Hz  on  the  Doppler  frequency  is 
used  to  discriminate  the  static  and  dynamic  part  of  the 
target.  The  spectrograms  of  the  resulting  radar  signals  are 
shown  in  Figs.  6a  and  6 b  for  the  rigid  body  and  the 
rotating  part,  respectively.  We  see  that  the  body  with 
nearly  constant  Doppler  and  the  rotating  part  with  fast 
changing  Doppler  are  separated. 

We  use  the  same  procedure  to  parameterise  and 
separate  the  radar  signals  from  range  cells  60  to  70. 
After  removing  the  rotating  part  interference,  the  final 
ISAR  image  is  shown  in  Fig.  7.  The  five  scatterers  of  the 
static  body  are  now  correctly  focused.  To  obtain  more 
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Fig.  5  Spectrogram  of  the  parameterised  radar  signal  using 
}00  chirplet  hoses 
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Fig.  6  Separated  body  and  rotating  part  signal 

a  Spectrogram  of  the  three  main  body  scallerers 
h  Spectrogram  of  the  rotating  part 
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Fig.  7  Reconstructed  ISAR  image  using  the  main  body  signals 
only 


time,  8 
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Fig.  8  Rotation  rate  estimation  using  autocorrelation 
a  Result  from  the  rotating  part  signal  after  separation 
b  Result  from  the  original  signal  before  separation 


information  about  the  rotation  motion,  an  autocorrelation 
analysis  of  the  separated  rotating  point  signal  is  shown  in 
Fig.  8a.  The  period  of  the  rotation  motion  is  determined  to 
be  0.15  s  from  this  Figure.  This  agrees  with  the  true 
rotation  rate  of  6.67  Hz.  For  comparison,  the  autocorrela¬ 
tion  of  the  raw  radar  signal  before  separation  is  shown  in 
Fig.  86.  It  is  difficult  to  detect  the  periodicity  from  the 
plot  as  the  rotating  scatterer  signal  is  heavily  contami¬ 
nated  by  the  large  body  return. 


4  Measurement  data  results 

The  algorithm  is  next  applied  to  two  sets  of  measurement 
data.  The  first  data  set  is  the  radar  data  collected  from  an  in¬ 
flight  aircraft  during  the  frontal  view  of  the  target.  The 
second  data  set  is  the  radar  data  collected  from  a  walking 
person.  In  both  cases,  the  goal  is  to  separate  the  main  body 
signal  from  the  rotating  part  return  for  better  target  feature 
extraction. 
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Fig.  9  !SAR  imaging  of  an  aircraft  during  frontal  view 


range 

Fig.  10  ISAR  image  of  the  aircraft  with  JEM  lines 


4. 1  Jet  engine  modulation  removal  from  an 
in-flight  aircraft 

The  geometry  of  the  problem  is  shown  in  Fig.  9.  The  radar 
collects  backscattering  data  from  an  in-flight  aircraft.  The 
resulting  ISAR  image  obtained  using  a  joint  time- 
frequency  based  motion  compensation  algorithm  [9]  is 
shown  in  Fig.  10.  We  observe  a  vertical  noisy  band  due  to 
the  rotating  engine  blades,  which  is  the  well  known  jet 
engine  modulation  (JEM)  phenomenon  [13].  The  geometry 
of  the  aircraft  body  is  obscured  due  to  the  presence  of  the 
JEM  lines. 

Simple  Doppler  gating  is  typically  used  to  alleviate  this 
problem.  The  result  in  Fig.  1 1  a  is  generated  by  putting  zeros 
in  cross-range  cells  1-32  and  62-128  in  the  image  area 
with  JEM  lines.  The  high  Doppler  frequency  components  in 
the  jet  engine  return  are  removed  in  this  manner.  However, 
we  see  that  in  areas  where  the  JEM  lines  overlap  with  the 
target  image,  this  technique  does  not  work  well,  as  it  cannot 
distinguish  the  aircraft  body  signal  from  the  JEM  signal 
with  low  Doppler  frequency. 

Using  the  chirplet-based  adaptive  signal  representation, 
we  first  parameterise  the  radar  signal.  Figure  11  h  is  the 
reconst  meted  ISAR  image  using  a  separation  criterion 
based  on  ft  only,  i.e.  we  remove  those  chirplet  bases  with 
large  ft  from  the  parameterised  signal.  It  is  much  better  than 
Fig.  1  \a  in  revealing  the  aircraft  body  feature.  This  confirms 
our  previous  observation  that  the  Doppler  rate  is  a  belter 
discriminator  than  the  Doppler  frequency  in  separating  the 
two  signals.  Finally,  we  use  both  ft  and  ft  to  separate  the 
two  signals.  The  aircraft  body  image  reconstructed  from 
chirplet  bases  with  both  small  ft  and  small  ft  is  shown  in 
Fig.  1 1  c.  We  see  an  even  better  representation  of  the  aircraft 


Fig.  11  Aircraft  body  and  JEM  line  separation 
a  Bodv  ISAR  image  with  Doppler  frequency  gating  only 
b  Bodv  ISAR  image  wilh  Doppler  rale  gating  only 

Body  ISAR  image  based  on  both  Doppler  frequency  and  Doppler  rate  parameters 

d  Separated  JEM  signal 
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body  feature  in  the  JEM  region.  The  JEM  signal  is  also 
displayed  in  Fig.  11  d.  The  signal  is  aliased  because  of  the 
low  PRF  of  the  radar  in  comparison  with  the  rotation  rate  ol 
the  en°ine  blades.  This  example  shows  that  this  algorithm 
works  despite  the  strong  Doppler  aliasing  of  the  rotating 
part  signal. 

4.2  Arm  swing  rate  estimation  from  a  walking 
person 

The  second  data  set  is  the  measured  radar  data  collected 
from  a  walking  person.  The  geometry  of  the  problem  is 
shown  in  Fig.  12.  Two  types  of  motions  are  involved:  the 
translation  motion  of  the  person’s  body  and  the  swinging 
motion  of  the  arms  (or  legs).  Figure  13  shows  the  range 
profiles  after  coarse  range  alignment  using  amplitude 
correlation.  Due  to  the  limited  range  resolution  relative  to 
the  tar°ct  size,  it  is  very  hard  to  discern  any  uselul  leatures 
about  either  the  body  or  the  arms  in  this  Figure.  Figure  14  is 
the  spectrogram  of  the  radar  signal  through  range  cell - 
Interesting  target  features  are  revealed  in  this  Figure.  The 
horizontal  Doppler  line  is  due  lo  the  body  motion  as  the 
person  walks  at  a  relatively  constant  speed  during  the  1-S  s 
dwell  interval.  The  sinusoidal-like  curve  shows  the  micro- 
Doppler  phenomenon  from  the  swinging  arm  motion.  The 
Doppler  spread  is  caused  by  the  varying  speed  ol  the  arm 
and  the  changing  angle  between  the  instantaneous  swinging 
motion  and  the  radar  incident  wave.  We  also  observe  the 
periodicity  of  the  arm  motion. 

To  separate  the  body  and  the  arm  returns,  the  adaptive 
chirplet  representation  is  applied.  After  the  parameterisa- 
tion.  we  again  separate  the  body  return  from  the  arm  return 
by  classifying  those  bases  with  large  Doppler  lrequency  Jk 
or  large  Doppler  rate  fa  as  contributions  from  the  arms.  The 
spectrograms  of  the  separated  body  and  arm  signals  are 
shown  in  Figs.  15 a  and  15/),  respectively.  The  main  features 
of  the  target  are  kept  after  the  separation,  indicating  good 
accuracy  of  the  parameterisation.  We  also  observe  a 
significant  denoising  efTect  from  the  parameterisation. 
This  is  because  noise  in  the  measured  data  does  not  have 
the  time -frequency  characteristics  of  a  chirplet  and  it  is  left 
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Fig.  13  Radar  data  after  range  compression 
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in  the  residual  signal  after  the  parameterisation.  The  arm¬ 
swing  period  can  be  easily  estimated  from  the  arm-only  data 
by  taking  the  autocorrelation  of  the  time  sequence.  The 
peaks  in  Fig.  16 a  correspond  to  the  period  of  the  signal, 
which  is  found  to  be  0.44  s.  Based  on  this  swing  rate  and  the 
speed  of  the  person  (2.3  m/s)  estimated  from  the  same  radar 
data,  the  stride  size  of  the  person  is  determined  to  be  about 
1.0  m.  For  comparison,  we  have  also  generated  the 
autocorrelation  of  the  original  data  without  the  joint 
lime- frequency  processing  (Fig.  16/?).  In  this  case, 
the  radar  return  from  the  arm  is  overshadowed  by  the  body 


dwell  time,  s 


Fig.  14  Spectrogram  of  the  radar  signal  containing  body  and 
ann  components 


Fig.  15  Separated  signal  components 
a  Body 

b  Swinging  arm 
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we  parameterise  and  separate  the  two  pans  using  the 
adaptive  signal  representation.  In  particular,  after  iormulat- 
ing  an  AM-FM  model  for  the  radar  signal  the  four- 
parameter  chi  rplet  basis  is  used  to  account  for  the  time  and 
frequency  localisation  of  the  signal.  After  the  parameterisa- 
tion.  the  separation  is  achieved  by  a  criterion  based  on  the 
extracted  Doppler  frequency  and  Doppler  rate  parameters. 
The  algorithm  has  been  successfully  tested  with  point- 
scatterer  simulations  and  applied  to  two  measurement  data 
sets.  In  the  aircraft  data,  we  are  able  to  reconstruct  a  better 
aircraft  body  image  after  the  separation.  In  the  walking 
person  data,  we  are  able  to  more  accurately  estimate  arm 
swing  rate.  The  results  demonstrate  the  potential  application 
of  this  algorithm  for  target  identification  using  1SAR  data 
from  non-rigid-body  targets. 
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Fig.  16  Arm-swing  rate  estimation  using  autocorrelation 

a  A  tier  signal  separation 
b  Before  signal  separation 

return  and  the  peaks  in  the  autocorrelation  function  are 
significantly  less  pronounced. 

Finally,  we  note  here  that  only  a  simple  exhaustive  search 
has  been  implemented  to  carry  out  the  chirplet  decompo¬ 
sition  in  our  examples.  In  the  walking  person  example,  the 
decomposition  of  a  signal  with  1024  data  samples  into 
50  chirplets  took  1050  s  using  MATLAB  on  a  personal 
computer  with  a  2.26 GHz  Pentium  4  CPU.  A  fast 
implementation  of  the  chirplet  decomposition  has  recently 
been  reported  in  [14]  and  should  speed  up  the  processing 
significantly. 

5  Conclusions 

In  this  paper,  a  chirplel-based  adaptive  signal  representation 
algorithm  has  been  applied  to  extract  features  from  1SAR 
data  of  a  target  with  a  rigid  main  body  and  a  rotating  part. 
Because  the  micro-Doppler  feature  of  the  rotating  part  is 
very  different  from  the  body  Doppler,  the  two  interfere  with 
each  other  if  processed  together.  To  overcome  this  problem. 
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A  formulation  for  ground  penetrating  radar  (GPR)  imaging  using  synthetic  aperture 
concept  is  introduced.  We  show  that  it  is  possible  to  form  a  3D  image  by  inverse  Fourier 
transforming  the  multi-frequency,  multi-spatial  scattered  field.  The  proposed  algorithm 
for  GPR  imaging  is  tested  with  measured  and  simulation  data. 
The  resulting  images  demonstrate  good  agreement  between  the  measured  and  simulated 
cases. 

Introduction:  The  imaging  of  buried  objects  or  inhomogeneities  underground  using  ground 
penetrating  radars  (GPR)  has  been  a  topic  of  interest  for  a  wide  variety  of  applications  ranging 
from  mine  detection  to  archeology.  Many  GPR  imaging  algorithms  have  been  proposed  in 
the  literature  [1-5],  Although  good  depth  resolution  can  usually  be  realized  in  GPR  images 
using  frequency  diversity,  good  resolution  in  the  cross-range  dimensions  is  much  harder  to 
achieve.  Capineri  et  al.  [3]  proposed  a  method  for  obtaining  good  resolution  in  GPR  images 
out  of  B-scan  data  by  applying  the  Hough  transformation  technique.  Moirow  and  Van 
Genderen  [4]  and  Van  Dongen  et  al.  [5]  applied  the  back  propagation  and  conjugate  gradient 
inversion  techniques  to  form  two-dimensional  (2D)  and  three-dimensional  (3D)  images  for  a 
borehole  radar.  However,  these  techniques  have  significant  computational  burden. 
Therefore,  there  is  a  need  for  obtaining  images  with  good  range  and  cross-range  resolution 
with  a  fast  algorithm. 

We  have  previously  developed  a  synthetic  aperture  algorithm  for  imaging  antenna- 
platform  interactions  based  on  multi-frequency,  multi-spatial  scattered  field  data  [6-8].  In  this 
paper,  we  extend  our  algorithm  to  generate  3D  GPR  images  of  scattered  data  from  buried 
objects  underground.  This  technique  is  based  on  the  approximate  Fourier  transform 


relationship  between  the  frequency-spatial  variables  and  the  distance-angle  information  of  the 
buried  scatterer.  The  algorithm  is  quite  attractive  since  it  forms  3D  images  by  using  a  fast 
Fourier  transform  (FFT)  followed  by  a  simple  transformation  from  the  distance-angle  domain 
to  the  image  domain.  It  is  computationally  fast.  Furthermore,  the  cross-range  resolution  can 
be  made  as  good  as  the  range  resolution  by  controlling  the  size  of  the  collection  aperture. 


SAR  Approach  for  GPR  Imaging:  Similar  to  the  antenna  synthetic  aperture  radar  algorithm  [6- 
8],  our  GPR  imaging  algorithm  is  based  on  collecting  the  multi-frequency  scattered  electric 
field  over  a  two-dimensional  spatial  grid  lying  on  top  of  the  ground  as  shown  in  Fig.  1.  We 
assume  the  target  point  P  is  located  at  an  unknown  location  fa,  yu  zf  We  also  assume  that 
the  frequency  bandwidth  is  small  compared  to  the  center  frequency  and  that  the  aperture 
dimensions  are  small  compared  to  Rn,  the  path  length  from  P  to  the  receiver.  Under  these 
assumptions,  the  scattered  electric  field  at  the  receiver  can  be  approximated  as  follows: 


Es{k,x',z')  =  Are 


-MRU+R2i) 


-  jkc  *'sin  or,-  .  p-jkc  -z'sin  or,  sin  fit  (1 ) 


where  A>  is  the  strength  of  the  scattered  field,  k  =  2nJe^U  is  the  wave  number  in  the  soil 
and  eris  the  relative  permittivity  of  the  ground.,  k  is  proportional  to  the  radar  frequency,  and 
kc  corresponds  to  the  wave  number  at  the  center  frequency.  By  taking  the  3D  inverse  Fourier 
transform  of  the  scattered  electric  field  with  respect  to  k,  x’  and  z\  it  is  possible  to  pinpoint 
the  total  travel  distance  and  the  angles  related  to  the  scatterer  location  as  follows: 

Es  (R,U,V)  =  Aj  ■  5{R  -  (R\i  +  %))  *  S(U  -  sin  at\ •  S(V  -  sin  at  sin  pj)  (2) 


Here,  we  introduce  three  new  variables  R=Ri+R2,  U=sina  and  F=sina»sin/?  for  simplicity. 
Once  an  image  in  the  (R,  U,  V)  domain  is  generated,  we  can  then  transform  it  from  the  (R,  U, 
V)  into  the  spatial  (x,  y,  z)  domain  by  using  the  trigonometric  relationship  between  the 


variables  (R,  a,  P)  and  (x,  y,  z).  The  transformations  from  (R,  JJ,  V)  to  (x,  y,  z)  is  unique  and 
correctly  maps  the  scatterer  location  [8].  However,  the  resultant  point  spread  response  in  the 
image  is  slightly  distorted  due  to  the  nonlinear  nature  of  the  transformation. 

Experimental  Results:  To  test  our  GPR  imaging  algorithm,  we  built  an  experimental  setup 
shown  in  Fig.  2.  In  this  setup,  a  wooden  pit  was  constructed  and  was  filled  with  play  sand. 
The  dielectric  constant  of  the  sand  was  measured  by  comparing  the  phase  delay  between  a 
pair  of  antennas  in  the  air  to  that  in  the  sand.  The  dielectric  constant  of  the  sand  was  found  to 
be  nearly  constant  at  2.26  for  the  frequency  range  from  5GHz  to  6GHz.  For  our  GPR 
experiment,  a  rectangular  copper  plate  whose  dimensions  are  46cm  in  the  x-direction  and 
30cm  in  the  z-direction  was  buried  at  46cm  below  the  sand  surface.  The  plate  is  located  at 
50cm  away  from  the  transmitter  along  the  x-axis.  The  S21  between  the  transmitter  and  the 
receiver  was  measured  using  an  HP8753C  network  analyzer.  As  the  transmitter  and  the 
receiver  antennas,  identical  coax-fed,  rectangular  waveguide  antennas  whose  dimensions  are 
3.81cm  and  1.91cm  were  used.  The  transmitter  antenna  was  assumed  to  be  placed  at  the 
origin  and  the  receiving  grid  was  assumed  to  be  centered  at  1  m  along  the  positive  x-direction. 
Both  antennas  were  horizontally  polarized  such  that  the  electric  field  was  parallel  to  the  metal 
plate.  The  scattered  field  was  collected  over  100  different  spatial  points.  The  size  of  the  10- 
by-10  receiving  grid  was  31.04cm  in  the  x-direction  and  14.83cm  in  the  z-direction.  For  all 
100  points,  the  signal  frequency  varied  from  4.9226  GHz  to  5.9352  GHz  over  25  evenly 
sampled  points.  The  measured  data  from  the  experimental  setup  in  Fig.2  were  collected  onto  a 
computer  and  processed.  After  applying  the  proposed  algorithm,  we  generated  a  3D  GPR 
image  of  the  region  below  the  surface.  Fig.3  (a)  shows  the  2D  projected  GPR  images  onto  the 
principal  Z-Y,  X-Y  and  X-Z  planes.  Overlaid  on  the  images  are  the  projected  outlines  of  the 
plate.  We  observe  two  main  hot  spots  in  the  image.  The  stronger  one  corresponds  to  the 


scattering  from  the  middle  of  the  plate  where  a  specular  point  exists.  The  weaker  one 
corresponds  to  the  diffraction  mechanism  from  the  front  edge  of  the  plate.  Both  image 
features  agree  well  with  the  geometrical  locations  of  the  plate.  In  addition,  resolutions  in  the 
cross  range  directions  (i.e.,  in  the  x-  and  z-directions)  are  nearly  the  same  as  the  resolution  in 
the  range  (y)  direction.  Note  that  the  spots  in  the  image  do  not  have  a  simple  point  spread 
form  and  they  are  somewhat  defocused.  This  is  due  the  non-linear  transformation  from  the  (R, 
a,  ft)  to  the  (x,  y,  z)  domain,  and  a  method  to  overcome  this  effect  has  been  discussed  in  [7, 8]. 
Nonetheless,  we  can  still  see  the  separation  between  the  two  points  on  the  plate,  which  are 
spaced  23  cm  apart  Therefore,  our  technique  is  able  to  achieve  good  resolutions  in  both  the 
range  and  the  cross-range  dimensions. 

Simulation  was  also  carried  out  using  a  physical  optics  calculation.  After  obtaining  the 
simulation  data  of  the  experimental  setup,  we  applied  the  same  imaging  algorithm  to  form  the 
simulated  GPR  image.  Fig.3  (b)  demonstrates  the  2D  projected  GPR  images  from  the 
simulation  data.  By  comparing  the  measured  GPR  images  to  the  simulated  ones,  we  see  good 
agreement  between  the  two.  Since  the  physical  theory  of  diffraction  contribution  was  absent 
in  the  simulation,  we  notice  that  the  edge  diffraction  contribution  in  the  simulation  is  weaker 
than  that  from  the  measured  image.  Finally,  data  were  also  collected  and  images  formed 
using  other  non-metallic  objects.  High-resolution  images  could  be  formed  consistently  using 
the  algorithm. 

Conclusion:  We  presented  a  Fourier  based  imaging  algorithm  for  ground  penetrating  radar 
based  on  the  synthetic  aperture  radar  concept.  The  algorithm  uses  the  phase  information  of  the 
scattered  field.  By  inverse  Fourier  transforming  the  scattered  field  data,  we  have  shown  that  it 
is  possible  to  form  high-resolution  3D  GPR  image  of  the  region  below  the  ground  surface.  To 
test  our  imaging  algorithm,  data  were  collected  from  a  buried  metallic  plate  using  an 


experimental  setup.  Our  imaging  algorithm  successfully  formed  a  3D  GPR  image  of  the  plate. 
The  measured  image  was  also  compared  to  that  formed  from  simulation  data  generated  using 
the  physical  optics  calculation.  Good  agreement  between  the  measured  and  simulated  images 
was  observed.  The  limitation  of  the  present  imaging  algorithm  is  that  it  assumes  the  soil 
medium  to  be  homogeneous  and  the  soil  property  is  known  a  priori. 
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